figure table

# Fabrication and characterization of optical nanofiber tips

Resmi M  School of Physics, University of Hyderabad, Hyderabad, Telangana,
India, 500046  Elaganuru Bashaiah  School of Physics, University of Hyderabad,
Hyderabad, Telangana, India, 500046  Ramachandrarao Yalla  School of Physics,
University of Hyderabad, Hyderabad, Telangana, India, 500046

###### Abstract

We demonstrate the fabrication of an optical nanofiber tip (ONFT) using a two-
step chemical etching technique. This technique employs 30% and 24%
hydrofluoric (HF) acid for the first and second steps, respectively. In the
first step, a silica single-mode fiber with clad diameter of 125 $\mu m$ and
core diameter of 10 $\mu m$ is immersed in the HF acid for 90 minutes. The
resultant fiber diameter is reduced to $\sim$45 $\mu m$. In the second step,
such a fiber is again immersed in the HF acid for 60-80 minutes. The etching
time is controlled to achieve the desired tip diameter of the ONFT. We
characterize fabricated ONFTs by measuring optical transmission and surface
morphology. The observed optical transmissions are more than 30% and tip
diameters are less than 500 $nm$. The measured results readily reveal the
merit of the employed technique as it is easy and cost-effective. Due to the
strong confinement of the electromagnetic field at the tip of ONFTs, these
structures lay versatile platforms with potential applications in sensing,
photonics, and quantum optics.

###### keywords:

Optical nanofiber tip, Chemical etching, Nanophotonics, Quantum optics

*Ramachandrarao Yalla<EMAIL_ADDRESS>
##  1 Introduction

A tapered optical fiber with a sub-wavelength-diameter tip is defined as an
optical nanofiber tip (ONFT). Due to high contrast between effective
refractive indices of its core (silica) and clad (air), the ONFT offers strong
field confinement [1, 2], leading to wide range of applications in fields of
photonics [3, 4], nanophotonics [5, 6], quantum optics [7, 8, 9, 10], and
sensing [11, 12, 13, 14]. Nowadays, the ONFT has become a promising candidate
due to its utilization in various research fields. However, the fabrication of
a high optical-quality ONFT is still a challenging task. For the fabrication
of the ONFT and optical nanofibers, various techniques have been proposed and
demonstrated [15, 16, 17, 18, 19, 20]. Examples would include flame-brushing
and pulling [15], CO2 laser-assisted heating and pulling [16], electric arc
heating and one-sided pulling [17], laser-heated pulling and bending [18],
mechanical polishing [19], and chemical etching [20]. Additionally, a
combination of electric arc heated micro-pulling with chemical etching [21]
and melt-stretching with chemical etching [22] have been demonstrated
experimentally.

In the flame brushing and pulling method [15], Lee et al. have demonstrated
the fabrication of a tapered fiber length of 217 $mm$, diameter of 820 $nm$,
and transmission of $\sim$90% (at a wavelength of 1550 $nm$). In the CO2
laser-assisted heating and pulling [16], Ward et al. have demonstrated the
fabrication of a fiber tapered diameter of $3$ \- $4$ $\mu m$ with the
transmission of $\sim$93% (at a wavelength of 980 $nm$). In the electric arc
heating and one-sided pulling technique [17], the fabrication of a fiber tip
with a diameter of 1.4 $\mu m$ was demonstrated. In the laser-heated pulling
and bending technique [18], the fabrication of a fiber tip with a short
tapering length of $\sim$300 $\mu m$ and a diameter of $\sim$50 $nm$ was
demonstrated. In the chemical etching technique using HF acid [20], the
fabrication of the ONFT with a diameter of 592 $nm$ was demonstrated.

Furthermore, a few modifications have been demonstrated in the chemical
etching technique. The ONFT has been fabricated using HF acid and protective
silicone oil on it [23, 24]. A syringe pump regulates the meniscus height on
the etchant’s surface, resulting in the desired cone angle. Dynamic selective
etching was used to demonstrate a parabolic shape optical nanofiber probe and
such a fiber tip shows improved transmission [25]. Translating the fiber
vertically during the etching process can control the shape of the fiber optic
tips [26]. Researchers have used dynamic etching to fabricate a variety of tip
shapes [27]. Short and long tapered ONFTs have been fabricated by translating
the fiber vertically upwards and downwards, respectively.

Among the above-mentioned techniques, chemical etching is economical as no
sophisticated instruments are required. It is easily adaptable and short
tapered ONFTs can be experimentally realized. Additionally, mechanical
stability of the ONFT is required for its application in different research
fields for robust usage. For example, chemically etched fiber tips have been
widely used for sensing [7, 28, 29, 30, 31, 32, 33] and quantum optics
applications [34, 35]. The ONFT with long tapering is mechanically unstable
compared to a short tapered ONFT. Hence, the chemical etching technique for
fabricating the short tapered ONFT with a sub-wavelength diameter tip draws
attention. Although various groups worldwide have demonstrated the ONFT
fabrication, the systematic fabrication and characterization of the ONFT have
not been reported yet.

In this paper, we demonstrate the fabrication of the ONFT using a two-step
chemical etching technique. This technique employs 30% and 24% HF acid for the
first and second steps, respectively. In the first step, a silica single-mode
fiber (SMF) with 125 $\mu m$ clad and 10 $\mu m$ core diameters is immersed in
the HF acid for 90 minutes. The resultant fiber diameter is reduced to
$\sim$45 $\mu m$. In the second step, such a fiber is again immersed in the HF
acid for 60 - 80 minutes. The measured optical transmissions are from 36% to
95%. The measured tip diameters of ONFTs are from 22 $nm$ to 820 $nm$.

##  2 Fabrication procedure

###  2.1 Pre-preparation process

Initially, we start with a 30 $cm$ long SMF (1550 BHP, Coherent) of single
mode cut-off wavelength at 1.3 $\mu m$. This fiber consists of a core, clad,
and outer plastic jacket with diameters 10 $\mu m$, 125 $\mu m$, and 250 $\mu
m$, respectively. Using a three-hole fiber stripping tool (FTS4, Thorlabs),
the outer jacket is carefully indented at one end, around 2 $cm$ from the
edge. Subsequently, the indented end is immersed in acetone for $5$ \- $10$
minutes to remove the outer plastic jacket without scratches on the fiber
surface. The resultant fiber diameter is reduced to 125 $\mu m$. The fiber is
meticulously cleaned with methanol to ensure it remains contaminant-free. The
stripped end of the fiber is cleaved perpendicularly using a fiber cleaver
(XL411, Thorlabs). The fiber is securely attached to a metal holder using a
UV-curable glue. The metal holder consists of a V-groove on one end so the
excess glue does not reach the cleaved fiber end. A small drop of the UV
curable glue is deposited at the V-groove location. Subsequently, the UV light
is directed to fix the fiber on the metal holder. Additionally, a scotch tape
is fixed on the UV glue location for more stability. The metal holder with the
SMF affixed on it, is placed on an xyz-translation stage (TS 65, Holmarc) to
achieve precise sample movement. Note that the fiber end is positioned
downwards for immersing into the HF meniscus.

###  2.2 Methodology for etching process

Figure 1: (a) A conceptual diagram for the fabrication of the optical
nanofiber tip. SMF, TS, MH, MT, and TH denote single-mode fiber, translation
stage, metal holder, microcentrifuge tube and tube holder, respectively. (b)
The picture of the experimental scenario.

A schematic of the experiment and the corresponding picture are shown in Figs.
1 (a) and (b), respectively. The whole experiment is set up in a fume
cupboard. A microcentrifuge tube of capacity 500 $\mu l$ is firmly fixed on a
tube holder. The original 40% HF acid (Avra, ASH2565) is diluted to 30% and
24% using de-ionized water. 30 $ml$ 40% HF and 10 $ml$ de-ionized water is
mixed to get 40 $ml$ of 30% HF acid. Also 24 $ml$ 40% HF and 16 $ml$ de-
ionized water is mixed to get 40 $ml$ of 24% HF acid. The translation stage is
placed near the microcentrifuge tube such that the SMF end is just above the
microcentrifuge tube. The microcentrifuge tube is filled with 30% and 24% HF
acid in the first and second steps, respectively. Then, HF acid forms a
meniscus on top of the microcentrifuge tube as clearly seen in Fig. 1 (b). The
length of the SMF dipped in the acid is precisely controlled by the
translation stage. In the first step, 2 $mm$ of the fiber end is immersed in
the acid for 90 minutes for the temperature ranging from 25-28∘C. The fiber is
carefully withdrawn from the acid at the end. Again, the fiber is immersed in
24% HF acid for 60/70/80 minutes. At the end, the fiber is withdrawn from the
acid, and any residual acid is rinsed off with de-ionized water to avoid any
further unintended etching of the fiber. The fabricated ONFT is preserved in a
clean environment to avoid any contamination that can result in degradation.

##  3 Characterization of the fabricated ONFT samples

We characterize optical and morphological properties of the fabricated ONFT
using laser and field emission scanning electron microscopy (FESEM) (Sigma 360
VP, Zeiss).

###  3.1 Optical

Figure 2: (a) A schematic of the experimental setup for optical transmission
measurements. LD, M, VND, FC, SMF, TS, MH, PD, C, ONFT, and DOPMC denote laser
diode, mirror, variable neutral density filter, fiber coupler, single mode
fiber, translation stage, metal holder, photodiode, connector, optical
nanofiber tip, and digital optical power meter console, respectively. (b) A
typical diffraction pattern as soon as the laser enters into guided modes of
the ONFT. (c) A typically measured input and output powers as a function of
time. (d) A histogram profile for the measured optical transmissions through
guided modes of ONFTs.

The experimental sketch is shown in Fig. 2 (a). The optical characterization
of the ONFT is performed by monitoring the power transmitted through its
guided modes. A laser diode (DL-G-5, Holmarc) at a wavelength of 532 $nm$ is
used. The laser light is directed through two mirrors to a variable neutral
density filter (HO-VND-N50, Holmarc) to vary the intensity of the laser light.
A SMF patch cable and a fiber coupler guide the laser into the ONFT. The SMF
patch cable has a single mode cut-off wavelength at 500-600 $nm$. Therefore it
filters all other modes and allows only the fundamental mode to the ONFT.
Single mode propogation is confirmed by observing the beam profile by
connecting the SMF patch cable to the laser, by connecting bare fiber to the
laser through the SMF patch cable before fabrication, and by connecting the
ONFT to the laser after fabrication. As shown in Fig. 2 (b), one can readily
observe a diffraction pattern. The optical power measurements are carried out
using a photodiode (S120VC, Thorlabs) connected to a digital optical power
meter console (PM100D, Thorlabs). The output power is measured for 15 minutes
by keeping the photodiode close to the tip of the ONFT. To measure the input
power, the ONFT is cleaved perpendicularly at a safe distance from the tip,
and the optical power at the cleaved end of the SMF is measured. Note that the
input power is measured as close to the ONFT to avoid coupling loss from the
SMF patch cable to the ONFT through the fiber coupler. The optical
transmission is derived as the ratio of the output power to the input power.
Figure 2 (c) shows a typically measured input and output powers as a function
of time. The black and red traces correspond to the output and input powers,
respectively. The measured transmission value is 75.0 ($\pm$1.5)%. We
performed such transmission measurements for several samples. A summary of
such measurements is plotted as a histogram and shown in Fig. 2 (d). The bin
width of the histogram is 5%, which is greater than the error in the
measurement. The histogram reflects that most ONFT samples exhibit optical
transmission $\geq$ 60%. The maximum/minimum measured optical transmission is
$\sim$95%/36%. The perpendicularly cleaved ONFT, after the output power
measurement, is used for the morphological characterization.

###  3.2 Morphological

Figure 3: Typical field emission scanning electron microscopy images of a few
optical fiber tips with different diameters. (a)-(b) and (c)-(f) are results
after the first and second steps, respectively. The measured diameters of (a),
(b), (c), (d), (e), and (f) are 45 ($\pm$1) $\mu m$, 47 ($\pm$2) $\mu m$, 820
($\pm$50) $nm$, 293 ($\pm$70) $nm$, 178 ($\pm$50) $nm$, and 22 ($\pm$7) $nm$,
respectively.

The ONFT is affixed to the metal holder using the UV curable glue as explained
in the section 2.1. Note that the ONFT is not free-hanging but rests on the
metal holder. To avoid charge-up, a mixture of gold and palladium is sputter
coated on samples for 3 minutes, resulting 2-3 $nm$ layer on the surface of
ONFTs. An accelerating voltage of 5 kV is applied for imaging. The tip of the
ONFT is magnified, and the diameter is measured. After the first step, typical
FESEM images of two samples are shown in Figs. 3 (a) and (b). The measured
diameters are 45 ($\pm$1) $\mu m$ and 47 ($\pm$2) $\mu m$ for (a) and (b),
respectively. After the second step, typical FESEM images of four ONFT samples
are shown in Figs. 3 (c)-(f). For (c), the measured ONFT diameter is 820
($\pm$50) $nm$ with etching timings of 90 and 70 minutes. The optical
transmission of 47.0 ($\pm$1.5)% is observed for it. For (d), the measured
ONFT diameter is 293 ($\pm$70) $nm$ with etching timings 60 and 30 minutes.
The observed optical transmission is 60 ($\pm$2)%. For (e), the measured ONFT
diameter is 178 ($\pm$50) $nm$ with etching timings 60 and 40 minutes. The
measured optical transmission is 81 ($\pm$1)%. For (f), the measured ONFT
diameter is 22 ($\pm$7) $nm$ with etching timings 90 and 80 minutes. The
optical transmission of 91 ($\pm$3)% is observed for it.

###  3.3 Data Analysis

Figure 4: Histogram profiles for measured diameters of optical nanofiber tips
with different etching timings of (a) 90 and 60 minutes (b) 90 and 70 minutes,
and (c) 90 and 80 minutes. Figure 5: (a) Typical tapering profiles as a
function of distance from the apex of ONFTs of tip diameters 183 $nm$ and 74
$nm$. Inset shows the tapering profile upto fiber diameter of 10 $\mu m$. (b)
Tapering angle as a function of distance from the apex of ONFTs with diameters
820 $nm$, 293 $nm$, and 178 $nm$, 183 $nm$ and 74 $nm$.

We measured ONFT diameters for several samples, and summaries for different
etching times are plotted as histograms as shown in Figs. 4 (a)-(c). Figure 4
(a) represents the histogram for ONFT diameters for the etching timings 90 and
60 minutes for the first and second steps, respectively. The minimum/maximum
measured ONFT diameter is 26/498 $nm$. Figure 4 (b) represents the histogram
for ONFT diameters for the etching timings 90 and 70 minutes for the first and
second steps, respectively. The minimum/maximum measured ONFT diameter is
27/820 $nm$. Figure 4 (c) represents the histogram for ONFT diameters for the
etching timings 90 and 80 minutes for the first and second steps,
respectively. The minimum/maximum ONFT diameter is 22/806 $nm$. Note that bin
widths of histograms are 50 $nm$, similar to measurement errors.

We also measured diameters of the ONFT at different locations in the tapering
region along the fiber axis. Typical tapering profiles as a function of
distance from the tip of the ONFT are shown in Fig. 5 (a). Note that tip
diameters are 183 $nm$ and 74 $nm$ and maximum diameter is 50.1 $\mu m$ for
both. The fiber diameter of 50 $\mu m$ is achieved at a distance of $\sim$0.5
$mm$ from the tip. We analysed diameter measurements to investigate tapering
angle. The tapering angle is calculated as the ratio of the difference in
diameters at two positions along the fiber axis and distance between them. The
summary of tapering angle as a function of the distance from the tip is shown
in Fig. 5 (b). The summary includes tapering angle profile for the three ONFTs
shown in Figs. 3 (c)-(e). and the two ONFTs mentioned in Fig. 5 (a). The
tapering angle varies from 0.08 $rad$ to 0.38 $rad$.

##  4 Discussions

As seen in Fig. 2 (b), the diffraction pattern is the first preview of the
fabricated ONFT. The measured result agrees with the previously reported
result in Ref. [25]. This may indicate the smooth tapering and the micro/nano
dimension of the ONFT. The measured optical transmission data from 36% to 95%
shown in Fig. 2 (d) indicates that ONFTs fabricated using this technique show
good optical quality. Most samples exhibit optical transmission of $\geq$ 60%,
confirming that the chemical etching employed here is an efficient technique
for fabricating optically good ONFTs. A summary of typical optical
transmissions with corresponding ONFT diameter can be found in Table. 1.
Diameter and optical transmission of previously reported results are also
shown for comparison. Note that the SMF patch cable used to guide the laser
light into the ONFT acts as a single mode filter. Therefore the measured
transmission is for the single mode. Hence higher transmission is observed. As
one can see in Table. 1, some samples show poor transmission. This may be due
to the multimode excitation in tapering region. Large variation of optical
transmissions might be due to the tapering profile and surface roughness.
Since the transmission is measured, it can be accounted in further experiments
for applications. We usually use optically good samples for further
experiments.

Table 1: Summary of typical ONFT diameters and corresponding optical transmission S. No. | Diameter of the ONFT ($nm$) | Optical transmission (%)
---|---|---
1 | 217 ($\pm$10) | 36 ($\pm$1)
2 | 204 ($\pm$10) | 46 ($\pm$2)
3 | 830 ($\pm$70) | 65 ($\pm$2.5)
4 | 274 ($\pm$20) | 90 ($\pm$1.5)
5 | 22 ($\pm$7) | 91 ($\pm$3)
6 | 820 ($\pm$50) | 47 ($\pm$1.5)
7 | 582 ($\pm$40) | 75 ($\pm$1.5)
8 | 127 ($\pm$15) | 95 ($\pm$1)
9 | 293 ($\pm$70) | 60 ($\pm$2)
10 | 178 ($\pm$50) | 81 ($\pm$1)
11 | 820 (Ref. [15]) | 90
12 | 3000 (Ref.[16]) | 93

After the first step, as seen in Figs. 3 (a) and (b), the measured diameters
are 45 $\mu m$ and 47 $\mu m$, respectively. The measured average etching rate
for the first step is 0.8 $\mu m/min$. As seen in Figs. 3 (c)-(f), the
measured diameters of ONFTs vary from 820 $nm$ to 22 $nm$. The measured
average etching rate for the second step is 0.55 - 0.75 $\mu m/min$. The
fluctuations in diameters may be due to the change in temperature and humidity
as these levels directly influence the etching rate. We recorded temperature
during the etching process and adjusted etching times acordingly. Although we
followed 90 minutes as the etching time for the first step for temperature in
the range of 25-28∘C, we reduced it to 60 minutes when the temperature was
$\geq$ 28∘C at the beginning of the first step. The etching time for the
second step was reduced to 30/40 minutes accordingly to achieve required fiber
diameter. Also, the SMF size was reduced to only a few micrometers in humid
condition and for temperature in the range of 23-25∘C. Note that the stability
of environmental parameters is required to ensure consistent results. As these
structures are so fragile, there is a chance that the tip may be broken while
withdrawing from the HF. Additionally, a small jerk also cause breaking of the
tip during experiment. At the end, the quality of the tip is asserted by
measuring the transmission and diffraction pattern for further experiments.

The fabrication of ONFTs with a wide range of diameters is confirmed using the
present technique. Note that ONFTs with minimum diameters of a few tens of
nanometers are experimentally realized. The trend of the histograms for three
particular etching times as shown in Figs. 4 (a)-(c) for the measured ONFT
diameters shows that the ONFT diameter does not vary as expected. We fabricate
ONFTs with diameters from tens to hundreds of nanometers for all the three
timings. The ONFT diameter does not decrease as the etching time is increased.
Though the ONFT diameters cannot be correlated to etching timings, the present
method is consistent to fabricate ONFTs of sizes less than 1 $\mu m$. At this
point of time, our lab is not equipped to investigate and precisely control
other parameters which affect the fabrication. In future we are planning to
control these parameters and to quantitatively study the effects. To have
precise control over the fabrication process, we are exploring heat and pull
technique [36, 37, 38] as another alternate approach. Consistent tip diameter
and transmission is promised in this technique.

The tapering profile shown in Fig. 5 (a), indicates that the fiber diameter
reaches 50 $\mu m$ at a distance of 0.5 $mm$ from the tip of the ONFT. We
dipped 2 $mm$ of the fiber in the HF for fabrication. This implies that the
actual fiber diameter of 125 $\mu m$ may be achieved within 2 $mm$.
Additionally, to improve the optical transmission with short tapering length,
we investigate different dipping lengths of the SMF in the HF acid. We dipped
the SMF into the HF acid up to 1.5 $mm$ and 1 $mm$. As the length dipped in
the HF gets reduced, shorter tapered ONFTs with good optical quality can be
realized. Such short-tapered ONFTs can be used in sensing and quantum optics
[39]. Tapering profiles for the two samples follow almost the same trend. It
indicates gradual tapering of the fiber, leading to high optical transmission.
The tapering angle shown in Fig. 5 (b), varies from 0.08 $rad$ to 0.38 $rad$.
It indicates that the tapering is not smooth compared to heat and pull
technique [36], in which tapering angle varies in $mrad$. Various samples do
not follow the same trend. Although it does not follow the trend as that of
heat and pull technique, chemically fabricated ONFTs are good enough to
perform experiments. Regarding the two-step etching, if 40% HF is used for the
entire etching process, the tapering will be very sharp and there will be
difficulty in controlling the size of the ONFT at nanoscales. If 24% HF is
used for the entire process, the etching process will take long time and the
environmental conditions will change considerably. Hence, a two-step etching
is preferred over a single step etching due to the control in etching rate.

Regarding the surface smoothness of the ONFT, as seen in Figs. 3, roughness
and crystalline structures are readily observed on the surface. Also, there is
a non-uniform etching throughout the fiber, especially in low temperatures.
Therefore, uncertainty in the ONFT diameter is inevitable in this technique.
To realize high-quality ONFTs, one can extend the fabrication of the ONFT
using the heat and pull technique [36, 37, 38]. In this technique, accurate
control over the ONFT diameter can be achieved. However, there is a trade-off
between the tip diameter and tapering length. Therefore, depending on the
application, one can adopt either chemical etching or heat and pull technique.

Although experimental ambiguities exist, we experimentally demonstrated the
fabrication of ONFTs with optical transmissions from 36% to 95% and tip
diameters ranging from 22 $nm$ to 820 $nm$. Note that always the tapered ONFT
is in submicron scale, which is suitable for applications. ONFTs sizes less
than 1 $\mu m$ is good enough for applications according to simulations [34,
40, 10, 41]. ONFTs with smaller diameters is efficient for sensing
applications [35, 12, 28, 29, 30, 13, 32]. Recently, our group reported
channeling of fluorescence photons from quantum dots into guided modes of the
ONFT which was fabricated using this technique [34]. The channeling was
confirmed by measuring the photon-counting and emission spectrum of
fluorescence photons. To calculate the channeling efficiency, optical
transmission is significant. Also measurement of tip diameter is essential to
determine whether the experiments are being carried out in the single mode
regime or multi mode regime. Hence ONFTs will find more applications in
quantum optics.

##  5 Conclusion

In summary, we demonstrated the development of an optical nanofiber tip (ONFT)
using a two-step chemical etching technique. We characterized fabricated ONFTs
by measuring optical transmission and surface morphology. The observed optical
transmissions are more than $30\%$ and tip diameters are less than $500$ $nm$.
Temperature, humidity and other environmental parameters have to be maintained
constant for consistent results. Consequently, tapering profile and surface
roughness can minimize variation in optical transmissions. The measured
results readily reveal the merit of the employed technique as easy and cost-
effective. Due to the strong confinement of the electromagnetic field at the
tip of ONFTs, these structures lay versatile platforms with potential
applications in sensing, photonics, and quantum optics.

### Disclosures

The authors declare no conflict of interests.

### Code, Data, and Materials Availability

Data may be obtained from authors upon reasonable request.

### Acknowledgments

This work was carried out in the wet chemistry lab under UGC-NRC in the School
of Physics, University of Hyderabad. The authors thank Shashank Suman and
Bratati Das for their helpful discussions. RM acknowledges the University
Grants Commission (UGC) for the financial support (Ref. No.:1412/CSIR-UGC NET
June 2019). RRY acknowledges the Science and Engineering Research Board (SERB)
for the Core Research Grant (CRG) (File No. CRG/2021/009185) and the Institute
of Eminence (IoE) grant at the University of Hyderabad, Ministry of Education
(MoE) (File No. RC2-21-019).

## References

  * [1] L. Tong, J. Lou, and E. Mazur, “Single-mode guiding properties of subwavelength-diameter silica and silicon wire waveguides,” Optics Express 12(6), 1025–1035 (2004).
  * [2] L. Tong and M. Sumetsky, Subwavelength and nanometer diameter optical fibers, Springer Science & Business Media (2011).
  * [3] Y. Xiong and F. Xu, “Multifunctional integration on optical fiber tips: challenges and opportunities,” Advanced Photonics 2(6), 064001–064001 (2020).
  * [4] X. Wu and L. Tong, “Optical microfibers and nanofibers,” Nanophotonics 2(5-6), 407–428 (2013).
  * [5] T. Tiecke, K. Nayak, J. D. Thompson, et al., “Efficient fiber-optical interface for nanophotonic devices,” Optica 2(2), 70–75 (2015).
  * [6] L. Tong, R. R. Gattass, J. B. Ashcom, et al., “Subwavelength-diameter silica wires for low-loss optical wave guiding,” Nature 426(6968), 816–819 (2003).
  * [7] J.-B. Decombe, S. Huant, and J. Fick, “Single and dual fiber nano-tip optical tweezers: trapping and analysis,” Optics express 21(25), 30521–30531 (2013).
  * [8] Y. Li, H. Xin, X. Xu, et al., “Fibre-optic trapping and manipulation at the nanoscale,” Advanced Materials Letters 9(8), 567–577 (2018).
  * [9] G. Brambilla, “Optical fibre nanowires and microwires: a review,” Journal of Optics 12(4), 043001 (2010).
  * [10] M. Resmi, E. Bashaiah, B. Das, et al., “Efficient fiber-coupled single photon source using an optical nanofiber tip,” in Women in Optics and Photonics in India 2022, 12638, 107–109, SPIE (2023).
  * [11] L. Tong, “Micro/nanofibre optical sensors: challenges and prospects,” Sensors 18(3), 903 (2018).
  * [12] S. Berneschi, A. Barucci, F. Baldini, et al., “Optical fibre micro/nano tips as fluorescence-based sensors and interrogation probes,” Optics 1(2), 213–242 (2020).
  * [13] J. S. Paiva, P. A. Jorge, C. C. Rosa, et al., “Optical fiber tips for biological applications: from light confinement, biosensing to bioparticles manipulation,” Biochimica et Biophysica Acta (BBA)-General Subjects 1862(5), 1209–1246 (2018).
  * [14] J. Lu and Z. Rosenzweig, “Nanoscale fluorescent sensors for intracellular analysis,” Fresenius’ Journal of Analytical Chemistry 366, 569–575 (2000).
  * [15] D. Lee, K. J. Lee, J.-H. Kim, et al., “Fabrication method for ultra-long optical micro/nano-fibers,” Current Applied Physics 19(12), 1334–1337 (2019).
  * [16] J. M. Ward, D. G. O’Shea, B. J. Shortt, et al., “Heat-and-pull rig for fiber taper fabrication,” Review of scientific instruments 77(8) (2006).
  * [17] W. Mukhtar, S. Shaari, and P. S. Menon, “Fabrication of optical fiber microprobe using electric arc heating and one-sided pulling technique,” in 2010 IEEE Student Conference on Research and Development (SCOReD), 104–106, IEEE (2010).
  * [18] A. Lazarev, N. Fang, Q. Luo, et al., “Formation of fine near-field scanning optical microscopy tips. part ii. by laser-heated pulling and bending,” Review of scientific instruments 74(8), 3684–3688 (2003).
  * [19] T. Grosjean, S. S. Saleh, M. A. Suarez, et al., “Fiber microaxicons fabricated by a polishing technique for the generation of bessel-like beams,” Applied optics 46(33), 8061–8067 (2007).
  * [20] H. J. Kbashi, “Fabrication of submicron-diameter and taper fibers using chemical etching,” Journal of Materials Science & Technology 28(4), 308–312 (2012).
  * [21] X. Huo, S. Pan, and S. Wu, “Fabrication of optical fiber probe nano-tips by heated micro-pulling combined with static chemical etching,” in 2006 1st IEEE International Conference on Nano/Micro Engineered and Molecular Systems, 254–257, IEEE (2006).
  * [22] H. Ren, C. Jiang, W. Hu, et al., “The preparation of optical fibre nanoprobe and its application in spectral detection,” Optics & Laser Technology 39(5), 1025–1029 (2007).
  * [23] H. Nikbakht, H. Latifi, M. Oraie, et al., “Fabrication of tapered tip fibers with a controllable cone angle using dynamical etching,” Journal of Lightwave Technology 33(23), 4707–4711 (2015).
  * [24] B. Puygranier and P. Dawson, “Chemical etching of optical fibre tips-experiment and model,” Ultramicroscopy 85(4), 235–248 (2000).
  * [25] W. Zhu, T. Shi, Z. Tang, et al., “Dynamic selective etching: a facile route to parabolic optical fiber nano-probe,” Optics Express 21(6), 6919–6927 (2013).
  * [26] L. Haber, R. Schaller, J. Johnson, et al., “Shape control of near-field probes using dynamic meniscus etching,” Journal of microscopy 214(1), 27–35 (2004).
  * [27] H. Muramatsu, K. Homma, N. Chiba, et al., “Dynamic etching method for fabricating a variety of tip shapes in the optical fibre probe of a scanning near-field optical microscope,” Journal of microscopy 194(2-3), 383–387 (1999).
  * [28] J. Ascorbe, J. Corres, I. Matias, et al., “High sensitivity humidity sensor based on cladding-etched optical fiber and lossy mode resonances,” Sensors and Actuators B: Chemical 233, 7–16 (2016).
  * [29] P. Zaca-Morán, J. Padilla-Martínez, J. Pérez-Corte, et al., “Etched optical fiber for measuring concentration and refractive index of sucrose solutions by evanescent waves,” Laser Physics 28(11), 116002 (2018).
  * [30] R. M. Chyad, M. Z. M. Jafri, and K. Ibrahim, “Fabrication nano fiber optic by chemical etching for sensing application,” Eng. Technol. 33, 994–1002 (2015).
  * [31] B. M. Cullum and T. Vo-Dinh, “The development of optical nanosensors for biological measurements,” Trends in biotechnology 18(9), 388–393 (2000).
  * [32] A. Pathak, V. Bhardwaj, R. Gangwar, et al., “Fabrication and characterization of tio2 coated cone shaped nano-fiber ph sensor,” Optics Communications 386, 43–48 (2017).
  * [33] S. K. Mondal, S. S. Pal, and P. Kapur, “Optical fiber nano-tip and 3d bottle beam as non-plasmonic optical tweezers,” Optics Express 20(15), 16180–16185 (2012).
  * [34] R. M, E. Bashaiah, and R. Yalla, “Channeling of fluorescence photons from quantum dots into guided modes of an optical nanofiber tip,” Journal of Optics 26, 065401 (2024).
  * [35] S. Suman, E. Bashaiah, R. M, et al., “In situ characterization of optical micro/nano fibers using scattering loss analysis,” Journal of Applied Physics 135, 123101 (2024).
  * [36] A. Stiebeiner, R. Garcia-Fernandez, and A. Rauschenbeutel, “Design and optimization of broadband tapered optical fibers with a nanofiber waist,” Optics express 18(22), 22677–22685 (2010).
  * [37] K. P. Nayak, M. Sadgrove, R. Yalla, et al., “Nanofiber quantum photonics,” Journal of Optics 20(7), 073001 (2018).
  * [38] J. Hoffman, S. Ravets, J. Grover, et al., “Ultrahigh transmission optical nanofibers,” AIP advances 4(6) (2014).
  * [39] K. G. Lee, X. W. Chen, H. Eghlidi, et al., “A planar dielectric antenna for directional single-photon emission and near-unity collection efficiency,” Nature Photonics 5(3), 166–169 (2011).
  * [40] S. Chonan, S. Kato, and T. Aoki, “Efficient single-mode photon-coupling device utilizing a nanofiber tip,” Scientific reports 4(1), 1–7 (2014).
  * [41] B. Das, M. Resmi, E. Bashaiah, et al., “Efficient single-mode coupling design using a silica/diamond nano-tip with a gold nanoparticle,” in Women in Optics and Photonics in India 2022, 12638, 125–127, SPIE (2023).

Ramachandrarao Yalla is an assistant professor at the University of Hyderabad,
India. He received his MSc in Physics from the University of Hyderabad, India,
in 2008 and his Ph.D. in Physics from the University of Electro-
Communications, Japan, in 2012. He is the author of more than 15 journal
papers, including Physical Review Letters and Applied Physics Letters. His
current research interests include quantum optics, nanophotonics, and
optoelectronic systems.

Resmi M is a PhD scholar at the University of Hyderabad, India. She received
her MSc degree in Physics from University of Madras, India in 2016. She is the
author of 5 articles including journal paper and conference paper.

Elaganuru Bashaiah is a PhD scholar at the University of Hyderabad, India. He
received his MSc degree in Physics from Vikrama Simhapuri University, India in
2016. He is the author of 5 articles including journal and conference papers.


# Enhancing SfePy with Isogeometric Analysis

Robert Cimrman∗† * Corresponding author<EMAIL_ADDRESS>New Technologies
Research Centre, University of West Bohemia, Plzeň, Czech RepublicCopyright ©
2014 Robert Cimrman. This is an open-access article distributed under the
terms of the Creative Commons Attribution License, which permits unrestricted
use, distribution, and reproduction in any medium, provided the original
author and source are credited. http://creativecommons.org/licenses/by/3.0/

###### Abstract

In the paper a recent enhancement to the open source package SfePy (Simple
Finite Elements in Python, http://sfepy.org) is introduced, namely the
addition of another numerical discretization scheme, the isogeometric
analysis, to the original implementation based on the nowadays standard and
well-established numerical solution technique, the finite element method. The
isogeometric removes the need of the solution domain approximation by a piece-
wise polygonal domain covered by the finite element mesh, and allows
approximation of unknown fields with a higher smoothness then the finite
element method, which can be advantageous in many applications. Basic
numerical examples illustrating the implementation and use of the isogeometric
analysis in SfePy are shown.

###### Index Terms:

partial differential equations, finite element method, isogeometric analysis,
SfePy

##  1 Introduction

Many problems in physics, biology, chemistry, geology and other scientific
disciplines can be described mathematically using a partial differential
equation (PDE) or a system of several PDEs. The PDEs are formulated in terms
of unknown field variables or fields, defined in some domain with a
sufficiently smooth boundary embedded in physical space.

SfePy (Simple Finite Elements in Python, http://sfepy.org) is a framework for
solving various kinds of problems (mechanics, physics, biology, …) described
by PDEs in two or three space dimensions. Because only the most basic PDEs on
simple domains (circle, square, etc.) can be solved analytically, a numerical
solution scheme is needed, involving, typically:

  * •

an approximation of the original domain by a polygonal domain;

  * •

an approximation of continuous fields by discrete fields defined by a finite
set of degrees of freedom (DOFs) and a (piece-wise) polynomial basis.

The above steps are called _discretization_ of the continuous problem. In the
following text two discretization schemes will be briefly outlined:

  * •

the finite element method [FEM] \- a long-established industry approved method
based on piece-wise polynomial approximation,

  * •

the isogeometric analysis [IGA] \- a quite recent generalization of FEM that
uses spline- or NURBS-based approximation.

SfePy, as its name suggests, has been based on FEM from its very beginning.
The IGA implementation has been added mainly due to the following reasons
(both will be addressed more in the text):

  * •

The IGA approximation can be globally smooth on a single patch geometry. The
continuity is determined by a few well defined parameters. This fact was the
main factor in deciding to implement IGA, because the smoothness is crucial in
one of our research applications (ab-initio electronic structure calculations
- work in progress). The high smoothness is paid for by the higher
computational complexity of the NURBS basis evaluation and higher fill-in of
the sparse matrix that a problem discretization leads to.

  * •

IGA can work directly with the geometric description of objects used in
geometric modeling and computer-aided design (CAD) systems, removing thus the
meshing step.

The paper is structured as follows. The geometric representation of objects is
outlined in Geometry Description using NURBS, because the terms defined there
are used in the IGA part of Outline of FEM and IGA. Then the particular
choices made in SfePy are presented in IGA Implementation in SfePy and
illustrated using examples of PDE solutions in Examples. All computations
below were done in SfePy, version 2014.3 - note that this paper is a short
description of the state and capabilities of the code as of this version. The
examples of numerical solutions have no particular scientific meaning or
importance besides being an illustration of the used methods.

##  2 Geometry Description using NURBS

First, let us briefly review the geometric representation of objects using
Bézier curves, B-splines and [NURBS] (Non-uniform rational B-spline) curves
and 2D (surface) or 3D (solid) bodies, to elucidate terminology used in
subsequent sections. Our IGA implementation is based on the explanation and
algorithms in [BE], thus below we follow its notation and definitions.

###  2.1 Bézier Curves

A Bézier curve is a parametric curve frequently used in computer graphics and
related fields. We define it here because its polynomial basis is used in the
code by means of the Bézier extraction technique, see [BE] and below.

A degree $p$ Bézier curve is defined by a linear combination of $p+1$
_Bernstein polynomial basis_ functions as

$C(\xi)=\sum_{a=1}^{p+1}\bm{P}_{a}B_{a,p}(\xi)=\bm{P}^{T}\bm{B}(\xi)\mbox{ for
}\xi\in[0,1]\;,$

where $\bm{P}=\\{\bm{P}_{a}\\}_{a=1}^{p+1}$ is the set of control points and
$\bm{B}(\xi)=\\{B_{a,p}(\xi)\\}_{a=1}^{p+1}$ is the set of Bernstein
polynomial basis functions. The Bernstein basis can be defined recursively for
$\xi\in[0,1]$ as $B_{a,p}(\xi)=(1-\xi)B_{a,p-1}(\xi)+\xi B_{a-1,p-1}(\xi)$,
$B_{1,0}(\xi)\equiv 1$, $B_{a,p}(\xi)\equiv 0$ if $a<1$ or $a>p+1$. An example
of Bézier curve is shown in Fig. 1.

Figure 1: A Bézier curve (blue) of degree three with four control points in
two space dimensions.

###  2.2 B-spline Curves

A B-spline is a generalization of the Bézier curve. B-splines (and their NURBS
generalization, see below) are used in computer graphics, geometry modeling
and related fields as well as the Bézier curves. In IGA B-spline basis
functions can be used for approximation of the unknown fields.

A univariate B-spline curve of degree $p$ is defined by a linear combination
of $n$ basis functions as

$T(\xi)=\sum_{A=1}^{n}\bm{P}_{A}N_{A,p}(\xi)=\bm{P}^{T}\bm{N}(\xi)\;,$

where $\bm{P}=\\{\bm{P}_{A}\\}_{A=1}^{n}$ is the set of control points. The
basis functions are defined by a _knot vector_ \- a set of non-decreasing
parametric coordinates $\Xi=\\{\xi_{1},\xi_{2},\dots,\xi_{n+p+1}\\}$, where
$\xi_{A}\in\mathbb{R}$ is the $A^{th}$ knot and $p$ is the polynomial degree
of the B-spline basis functions. Then for $p=0$

$\displaystyle N_{A,0}(\xi)$ $\displaystyle=$ $\displaystyle 1\mbox{ for
}\xi_{A}\leq\xi<\xi_{A+1}\;,$ $\displaystyle=$ $\displaystyle 0\mbox{
otherwise.}$

For $p>0$ the basis functions are defined by the Cox-de Boor recursion formula

$N_{A,p}(\xi)=\frac{\xi-\xi_{A}}{\xi_{A+p}-\xi_{A}}N_{A,p-1}(\xi)+\frac{\xi_{A+p+1}-\xi}{\xi_{A+p+1}-\xi_{A+1}}N_{A+1,p-1}(\xi)\;.$

Note that it is possible to insert knots into a knot vector without changing
the geometric or parametric properties of the curve by computing the new set
of control points in a particular way, see e.g. [BE].

A B-spline curve with a knot vector with no internal knots, i.e. of the form

$\Xi=\\{\underbrace{0,\dots,0}_{p+1},\underbrace{1,\dots,1}_{p+1}\\}\;,$

corresponds to a Bézier curve of degree $p$ with the same control points.

###  2.3 NURBS Curves

B-splines can be used to approximately describe almost any geometry. Their
main drawback is the fact, that a circular or spherical segment cannot be
described exactly. This problem was eliminated by the introduction of NURBS in
geometry modelling.

A NURBS (Non-uniform rational B-spline) of degree $p$ is defined by a linear
combination of $n$ rational basis functions as

$T(\xi)=\sum_{A=1}^{n}\bm{P}_{A}R_{A,p}(\xi)=\bm{P}^{T}\bm{R}(\xi)\;,$

where $\bm{P}=\\{\bm{P}_{A}\\}_{A=1}^{n}$ is the set of control points and
$\bm{R}(\xi)=\\{R_{A,p}(\xi)\\}_{A=1}^{p+1}$ is the set of rational basis
functions. The rational basis functions are defined using the B-spline basis
functions as

$R_{A,p}(\xi)=\frac{w_{A}{N_{A,p}(\xi)}}{W(\xi)}\;,\quad
W(\xi)=\sum_{B=1}^{n}w_{B}N_{B,p}(\xi)\;,$

where $w_{i}$ is the weight corresponding to the $i^{th}$ basis function and
$W$ is the weight function.

Note that a NURBS curve in $\mathbb{R}^{n}$ is equal to a B-spline curve in
$\mathbb{R}^{n+1}$:

$T(\xi)=\sum_{A=1}^{n}\bar{\bm{P}_{A}}N_{A,p}(\xi)\;,\quad\bar{\bm{P}_{A}}=\\{w_{A}\bm{P}_{A},w_{A}\\}^{T}\;.$

This means that all algorithms that work for B-splines work also for NURBS.

####  2.3.1 NURBS Surfaces and Solids

A surface is obtained by the tensor product of two NURBS curves. The knot
vector is defined for each axial direction and there are $n\times m$ control
points for $n$ basis functions in the first axis and $m$ basis functions in
the second one.

Analogically, a solid is given by tensor product of three NURBS curves.

####  2.3.2 NURBS Patches

Complex geometries cannot be described by a single NURBS outlined above, often
called _NURBS patch_ \- many such patches might be needed, and special care
must be taken to ensure required continuity along patch boundaries and to
avoid holes. A single patch geometry will be used in the following text, see
Fig. 2.

##  3 Outline of FEM and IGA

The two discretization methods will be illustrated on a very simple PDE - the
Laplace equation - in a plane (2D) domain. The Laplace equation describes
diffusion and can be used to determine, for example, temperature or electrical
potential distribution in the domain. We will use the "temperature"
terminology and the notation from Table I.

symbol |  meaning
---|---
$\Omega$ |  solution domain
$\Omega_{h}$ |  discretized solution domain
$\Gamma_{D}$, $\Gamma_{N}$ |  subdomains representing parts of the domain surface for applying Dirichlet and Neumann boundary conditions
$\underline{n}$ |  unit outward normal
$\nabla\equiv[\frac{\partial}{\partial x_{1}},\frac{\partial}{\partial x_{2}}]^{T}$ |  gradient operator
$\nabla\cdot$ |  divergence operator
$\Delta\equiv\nabla\cdot\nabla$ |  Laplace operator
$C^{1}$ |  space of functions with continuous first derivatives
$H^{1}$ |  space of functions with integrable values and first derivatives
$H^{1}_{0}$ |  space of functions from $H^{1}$ that are zero on $\Gamma_{D}$
TABLE II: Notation.

The problem is as follows: Find temperature $T$ such that:

$\displaystyle\Delta T$ $\displaystyle=$ $\displaystyle 0\mbox{ in }\Omega\;,$
(1) $\displaystyle T$ $\displaystyle=$ $\displaystyle\bar{T}\mbox{ on
}\Gamma_{D}\;,$ (2) $\displaystyle\nabla T\cdot\underline{n}$ $\displaystyle=$
$\displaystyle 0\mbox{ on }\Gamma_{N}\;,$ (3)

where the second equation is the Dirichlet (or essential) boundary condition
and the third equation is the Neumann (or natural) boundary condition that
corresponds to a flux through the boundary.

The operator $\Delta$ has second derivatives - that means that the solution
$T$ needs to have continuous first derivatives, or, it has to be from $C^{1}$
function space - this is often not possible in examples from practice.
Instead, a _weak solution_ is sought that satisfies: Find $T\in H^{1}(\Omega)$

$\displaystyle\int_{\Omega}\nabla s\cdot\nabla
T-\underbrace{\int_{\Gamma_{N}}s\ \nabla T\cdot\underline{n}}_{\equiv 0}$
$\displaystyle=$ $\displaystyle 0\;,\quad\forall s\in H^{1}_{0}(\Omega)\;,$
(4) $\displaystyle T$ $\displaystyle=$ $\displaystyle\bar{T}\quad\mbox{ on
}\Gamma_{D}\;,$ (5)

where the natural boundary condition appears naturally in the equation (hence
its name). The above system can be derived by multiplying the original
equation by a test function $s\in H^{1}_{0}(\Omega)$, integrating over the
whole domain and then integrating by parts.

Both FEM and IGA now replace the infinite function space $H^{1}(\Omega)$ by a
finite subspace with a basis with a small support on a discretized domain
$\Omega_{h}$, see below for particular basis choices. Then
$T(\underline{x})\approx\sum_{k=1}^{N}T_{k}\phi_{k}(\underline{x})$, where
$T_{k}$ are the DOFs and $\phi_{k}$ are the base functions. Similarly,
$s(\underline{x})\approx\sum_{k=1}^{N}s_{k}\phi_{k}(\underline{x})$.
Substituting those into (4) we obtain

$\displaystyle\int_{\Omega_{h}}\left(\sum_{j=1}^{N}s_{j}\nabla\phi_{j}\cdot\sum_{k=1}^{N}\nabla\phi_{k}T_{k}\right)=0\;.$

This has to hold for any $s$, so we can choose $s=\phi_{j}$ for $j=1,\dots,N$.
It is also possible to switch the sum with the integral and put the constants
$T_{k}$ out of the integral, to obtain the discrete system:

$\displaystyle\sum_{k=1}^{N}\int_{\Omega_{h}}\left(\nabla\phi_{j}\cdot\nabla\phi_{k}\right)T_{k}=0\;.$
(6)

In compact matrix notation we can write $\bm{K}\bm{T}=\bm{0}$, where the
matrix $\bm{K}$ has components
$K_{ij}=\int_{\Omega_{h}}\nabla\phi_{i}\cdot\nabla\phi_{j}$ and $\bm{T}$ is
the vector of $T_{k}$. The Dirichlet boundary conditions are satisfied by
setting the $T_{k}$ on the boundary $\Gamma_{D}$ to appropriate values.

Both methods make use of the small support and evaluate (6) as a sum over
small "elements" to obtain local matrices or vectors that are then assembled
into a global system - system of linear algebraic equations in our case.

The particulars of domain geometry description and basis choice will now be
outlined. For both methods, we will use the domain shown in Figure 2. Its
geometry is described by NURBS, see Geometry Description using NURBS.

Figure 2: The domain with a curved boundary described by NURBS. The internal
lines correspond to several selected iso-lines given by the parametrization of
the 2D tensor-product NURBS patch.

###  3.1 FEM

In this method a continuous solution domain is approximated by a polygonal
domain - _FE mesh_ \- composed of small basic subdomains with a simple
geometric shape (e.g. triangles or quadrilaterals in 2D, tetrahedrons or
hexahedrons in 3D) - the elements. The continuous fields of the PDEs are
approximated by polynomials defined on the individual elements. This
approximation is (usually) continuous over the whole domain, but its
derivatives are only piece-wise continuous.

First we need to make a FE mesh from the NURBS description, usual in CAD
systems. While it is easy for our domain, it is a difficult task in general,
especially in 3D space. Here a cheat has been used and the mesh depicted in
Figure 3 was generated from the NURBS description using the IGA techniques
described below. Quite a fine mesh had to be used to capture the curved
boundaries.

Figure 3: The FE-discretized domain covered by quadrilateral elements, forming
the FE mesh.

Having the geometry discretized, a suitable approximation of the fields has to
be devised. In (classical1) FEM, the base functions with small support are
polynomials, see Figure 4 for an illustration in 1D. A $k$-th base function is
nonzero only in elements that share the DOF $T_{k}$ and it is a continuous
polynomial over each
element.\raisebox{10.00002pt}{\hypertarget{id8}{}}\hyperlink{id7}{1}\raisebox{10.00002pt}{\hypertarget{id8}{}}\hyperlink{id7}{1}footnotetext:
See the Wikipedia page for a basic overview of FEM and its many variations:
http://en.wikipedia.org/wiki/Finite_element_method.

Figure 4: The 1D FE basis on three line elements with black thick line an
interpolated function resulting from the same DOF vector for each row: top:
linear, bottom: quadratic, left: Lagrange, right: Lobatto. Each basis function
has a single color.

The thick black lines in Figure 4 result from interpolation of the DOF vector
generated by $\sin(\frac{\pi}{2}\frac{x}{3})$ evaluated in points of maximum
of each basis function. The left column of the figure shows the Lagrange
polynomial basis, which is interpolatory, i.e., a DOF value is equal to the
approximated function value in the point, called _node_ , where the basis is
equal to 1. The right column of the figure shows the Lobatto polynomial basis,
that is not interpolatory for DOFs belonging to basis functions with order
greater than 1 - that is why the bottom right interpolated function differs
from the other cases. This complicates several things (e.g. setting of
Dirichlet boundary conditions - a projection is needed), but the hierarchical
nature of the basis, i.e. increasing approximation order means adding new
basis functions without modifying the existing ones, has also advantages, for
example better condition number of the matrix for higher order approximations.

The basis functions are usually defined in a reference element, and are then
mapped to the physical mesh elements by an (affine) transformation. For our
mesh we will use bi-quadratic polynomials over the reference quadrilateral - a
quadratic function along each axis direction, such as those in the bottom row
of Figure 4.

Several families of the element basis functions exist. In SfePy, Lagrange
basis and Lobatto (hierarchical) basis can be used on quadrilaterals, see
Figure 5.

Figure 5: Bi-quadratic basis functions on the reference quadrilateral: left:
Lagrange right: Lobatto.

###  3.2 IGA

In IGA, the CAD geometrical description in terms of NURBS patches is used
directly for the approximation of the unknown fields, without the intermediate
FE mesh - the meshing step is removed, which is one of its principal
advantages. As described in Geometry Description using NURBS, a D-dimensional
geometric domain is defined by

$\underline{x}(\underline{\xi})=\sum_{A=1}^{n}\bm{P}_{A}R_{A,p}(\underline{\xi})=\bm{P}^{T}\bm{R}(\underline{\xi})\;,$

where $\underline{\xi}=\\{\xi_{1},\dots,\xi_{D}\\}$ are the parametric
coordinates, and $\bm{P}=\\{\bm{P}_{A}\\}_{A=1}^{n}$ is the set of control
points. The same NURBS basis is used also for the approximation of PDE
solutions. For our temperature problem we have

$T(\underline{\xi})=\sum_{A=1}^{n}T_{A}R_{A,p}(\underline{\xi})\;,\quad
s(\underline{\xi})=\sum_{A=1}^{n}s_{A}R_{A,p}(\underline{\xi})\;,$

where $T_{A}$ are the unknown DOFs - coefficients of the basis in the linear
combination, and $s_{A}$ are the test function DOFs.

Figure 6: From left to right: parametric mesh (tensor product of knot
vectors), control mesh, Bézier mesh.

Our domain in Figure 2 can be exactly described by a single NURBS patch.
Several auxiliary grids (called "meshes" as well, but do not mistake with the
FE mesh) can be drawn for the patch, see Figure 6. The parametric mesh is
simply the tensor product of the knot vectors defining the parametrization -
the lines correspond to the knot vector values. In our case there are four
unique knot values in the first parametric axis and five in the second axis.
The control mesh has vertices given by the NURBS patch control points and
connectivity corresponding to the tensor product nature of the patch. The
Bézier mesh will be described below. The thin blue lines are iso-lines of the
NURBS parametrization, as in Figure 2.

On a single patch, such as our whole domain, the NURBS basis can be
arbitrarily smooth - this is another compelling feature not easily obtained by
FEM. The basis functions $R_{A,p}$, $A=1,\dots,n$ on the patch are uniquely
determined by the knot vector for each axis, and cover the whole patch, see
Figure 7.

Figure 7: The degree 2 NURBS basis functions on the single patch domain.

##  4 IGA Implementation in SfePy

Our implementation uses a variant of IGA based on Bézier extraction operators
[BE] that is suitable for inclusion into existing FE codes. The code itself
does not see the NURBS description at all. The NURBS description can be
prepared, for example, using igakit package, a part of [PetIGA].

The Bézier extraction is illustrated in Figure 8. It is based on the
observation that repeating a knot in the knot vector decreases continuity of
the basis in that knot by one. This can be done in such a way that the overall
shape remains the same, but the "elements" appear naturally as given by non-
zero knot spans. The final basis restricted to each of the elements is formed
by the Bernstein polynomials $\bm{B}$.

In [BE] algorithms are developed that allow computing _Bézier extraction
operator_ $\bm{C}$ for each such element such that the original (smooth) NURBS
basis function $\bm{R}$ can be recovered from the local Bernstein basis
$\bm{B}$ using $\bm{R}=\bm{C}\bm{B}$. The Bézier extraction also allows
construction of the Bézier mesh, see Figure 6, right. The code then loops over
the Bézier elements and assembles local contributions in the usual FE sense.

Figure 8: From left to right: NURBS basis of degree 2 that describes the
second axis of the parametric mesh, corresponding Bernstein basis with Bézier
elements delineated by vertical lines.

In SfePy, various subdomains can be defined using _regions_ , see [SfePy]. For
example, below we use the following region definition to specify an internal
subdomain:

>
>     ’vertices in (x > 1.5) & (y < 1.5)’
>

To make this work with IGA, where no real mesh exists, a _topological Bézier
mesh_ is constructed, using only the corner vertices of the Bézier mesh
elements, because those are interpolatory, i.e., they are in the domain or on
its boundary, see Figures 6, 8 right.

The regions serve both to specify integration domains of the terms that make
up the equations and to define the parts of boundary, where boundary
conditions are to be applied. SfePy supports setting the Dirichlet boundary
conditions by user-defined functions of space (and time). To make this feature
work with IGA, a projection of the boundary condition functions to the space
spanned by the appropriate boundary basis functions was implemented.

###  4.1 Notes on Code Organization

Although the Bézier extraction technique shields the IGA-specific code from
the rest of the FEM package, the implementation was not trivial. Similar to
the Lobatto FE basis, the DOFs corresponding to the NURBS basis are not equal
to function values with the exception of the patch corners. Moreover, the IGA
fields do not work with meshes at all - they need the NURBS description of the
domain together with the Bézier extraction operators and the topological
Bézier mesh. So the original sfepy.fem sub-package was renamed and split into:

  * •

sfepy.discrete for the general classes independent of the particular
discretization technique (for example variables, equations, boundary
conditions, materials, quadratures, etc.);

  * •

sfepy.discrete.fem for the FEM-specific code;

  * •

sfepy.discrete.iga for the IGA-specific code;

  * •

sfepy.discrete.common for common functionality shared by some classes in
sfepy.discrete.fem and sfepy.discrete.iga.

In this way, circular import dependencies were minimized.

###  4.2 Using IGA

As described in [SfePy], problems can be described either using problem
description files - Python modules containing definitions of the various
components (mesh, regions, fields, equations, …) using basic data types such
as dict and tuple, or using the sfepy package classes directly interactively
or in a script. The former way is more basic and will be used in the
following.

In a FEM computation, a mesh has to be defined using:

filename_mesh = ’fe_domain.mesh’In an IGA computation, a NURBS domain has to
be defined instead:

filename_domain = ’ig_domain.iga’where the ’.iga’ suffix is used for a custom
HDF5 file that can be prepared by functions in sfepy.discrete.iga.

A scalar real FE field with the approximation order 2 called ’temperature’ can
be defined by:

# Lagrange basis is the default.fields = { ’temperature’ : (’real’, 1,
’Omega’, 2),}# Lobatto basis.fields = { ’temperature’ : (’real’, 1, ’Omega’,
2, ’H1’, ’lobatto’),}An analogical IGA field can be defined by:

fields = { ’temperature’ : (’real’, 1, ’Omega’, None, ’H1’, ’iga’),}Here the
approximation order is None, as it is given by the ’.iga’ domain file.

The above are the only changes required to use IGA - everything else remains
the same as in FEM calculations. The scalar and vector volume terms (weak
forms, linear or nonlinear) listed at http://sfepy.org/doc-
devel/terms_overview.html can be used without modification.

###  4.3 Limitations

There are currently several limitations that will be addressed in future:

  * •

projections of functions into the NURBS basis;

  * •

support for surface integrals;

  * •

linearization of results for post-processing;

    * –

currently the fields on a tensor-product patch are sampled by fixed parameter
vectors and a corresponding FE-mesh is generated;

  * •

all variables have to have the same approximation order, as the basis is given
by the domain file;

  * •

the domain is a single NURBS patch only.

##  5 Examples

Numerical examples illustrating the IGA calculations are presented below. The
corresponding problem description files for SfePy (version 2014.3) can be
downloaded from https://github.com/sfepy/euroscipy2014-iga-examples, revision
28b1fe9bff7043da6fd159c20b9f244337a17e82, or from
http://dx.doi.org/10.5281/zenodo.12257.

###  5.1 Temperature Distribution

The 2D domain depicted in Figure 2 is used in this example. The temperature
distribution is given by the solution of the Laplace equation (4) with a
particular set of Dirichlet boundary conditions on $\Gamma_{D}$. The region
$\Gamma_{D}$ consisted of two parts $\Gamma_{1}$, $\Gamma_{2}$ of the domain
boundary on the opposite edges of the patch, see Figure 9 \- the temperature
was fixed to 0.5 on $\Gamma_{1}$ and to -0.5 on $\Gamma_{2}$, as can be seen
in Figure 10. As mentioned in Limitations, the resulting field $T$ was sampled
by fixed uniform parameter vectors along each axis, and the corresponding
output FE mesh was generated. The mesh was saved in the VTK format and the
results visualized using SfePy’s postproc.py script based on Mayavi. The
generated mesh can be seen as the undeformed wire-frame.

Figure 9: The regions defined on the domain shown on the topological Bézier
mesh by red color. From left: $\Gamma_{1}$, $\Gamma_{2}$, $\Omega_{0}$

For comparison with a FEM solution, see Figure 11. The FEM problem had 1363
DOFs in the linear system, while the IGA problem only 20. The mesh depicted in
Figure 3 was used for the FEM computation.

Figure 10: A solution of the 2D Laplace equation.

Figure 11: A solution of the 2D Laplace equation by FEM.

Next we added a negative source term to the Laplace equation in region
$\Omega_{0}$ (see Figure 9 right) to obtain the Poisson equation:

$\displaystyle\int_{\Omega}\nabla s\cdot\nabla T$ $\displaystyle=$
$\displaystyle\int_{\Omega_{0}}-2s\;,\quad\forall s\in H^{1}_{0}(\Omega)\;,$
(7) $\displaystyle T$ $\displaystyle=$ $\displaystyle\bar{T}\quad\mbox{ on
}\Gamma_{D}\;,$ (8)

The corresponding solution can be seen in Figure 12. The boundary conditions
stayed the same as in the previous case.

Figure 12: A solution of the 2D Poisson equation with volume source in a
subdomain.

The complete problem description file for computing (7) is shown below. See
[SfePy] or http://sfepy.org for explanation.

filename_domain = ’ig_domain.iga’regions = { ’Omega’ : ’all’, ’Omega_0’ :
’vertices in (x > 1.5) & (y < 1.5)’, ’Gamma1’ : (’vertices of set xi10’,
’facet’), ’Gamma2’ : (’vertices of set xi11’, ’facet’),}fields = {
’temperature’ : (’real’, 1, ’Omega’, None, ’H1’, ’iga’),}variables = { ’T’ :
(’unknown field’, ’temperature’, 0), ’s’ : (’test field’, ’temperature’,
’T’),}ebcs = { ’T1’ : (’Gamma1’, {’T.0’ : 0.5}), ’T2’ : (’Gamma2’, {’T.0’ :
-0.5}),}materials = { ’m’ : ({’f’ : -2.0},),}integrals = { ’i’ : 3,}equations
= { ’Temperature’ : """dw_laplace.i.Omega(s, T) = dw_volume_lvf.i.Omega_0(m.f,
s)"""}solvers = { ’ls’ : (’ls.scipy_direct’, {}), ’newton’ : (’nls.newton’, {
’i_max’ : 1, ’eps_a’ : 1e-10, }),}

###  5.2 Elastic Deformation

This example illustrates a calculation with a vector variable, the
displacement field $\underline{u}$, given by deformation of a 3D elastic body.
The weak form of the problem is: Find $\underline{u}\in[H^{1}(\Omega)]^{3}$
such that:

$\displaystyle\int_{\Omega}D_{ijkl}\
e_{ij}(\underline{v})e_{kl}(\underline{u})$ $\displaystyle=$ $\displaystyle
0\;,\quad\forall\underline{v}\in[H^{1}_{0}(\Omega)]^{3}\;,$
$\displaystyle\underline{u}$ $\displaystyle=$
$\displaystyle\bar{\underline{u}}\quad\mbox{ on }\Gamma_{D}\;,$

where
$D_{ijkl}=\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+\lambda\delta_{ij}\delta_{kl}$
is the isotropic stiffness tensor given in terms of Lamé’s coefficients
$\lambda$, $\mu$ and $e_{ij}(\underline{u})=\frac{1}{2}(\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})$ is the Cauchy,
or small strain, deformation tensor. The equation expresses the internal and
external (zero here) force balance, where the internal forces are described by
the Cauchy stress tensor $\sigma_{ij}(\underline{u})=D_{ijkl}\
e_{kl}(\underline{u})$.

The 3D domain $\Omega$ was simply obtained by extrusion of the 2D domain of
the previous example, and again $\Gamma_{D}$ consisted of two parts
$\Gamma_{1}$, $\Gamma_{2}$. The body was clamped on $\Gamma_{1}$:
$\underline{u}=0$ and displaced on $\Gamma_{2}$: $u_{1}=0.01$,
$u_{2}(\underline{x})=-0.02x_{2}$ and
$u_{3}(\underline{x})=-0.02+(0.15*(x_{1}-1)^{2})$, for
$\underline{x}\in\Gamma_{2}$. Note that the Dirichlet boundary conditions on
$\Gamma_{2}$ depend on the position $\underline{x}$. The corresponding
solution can be seen in Figure 13.

Figure 13: A solution of the 3D linear elasticity equation. The undeformed
domain is shown as a wireframe, 10x magnified deformation.

##  6 Conclusion

Two numerical techniques for discretization of partial differential equations
were briefly outlined and compared, namely the well-established and proven
finite element method and its much more recent generalization, the
isogeometric analysis, on the background given by the open source finite
element package SfePy, that has been recently enhanced with the isogeometric
analysis functionality.

The Bézier extraction operators technique, that was used for a relatively
seamless integration into the existing finite element package, was mentioned,
as well as some of the difficulties "on the road" and limitations of the
current version.

Numerical examples - a scalar diffusion problem in 2D and a vector elastic
body deformation problem in 2D were shown.

###  6.1 Support

Work on SfePy is partially supported by the Grant Agency of the Czech
Republic, project P108/11/0853.

## References

  * [FEM] Thomas J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, 2000.
  * [IGA] J. Austin Cottrell, Thomas J.R. Hughes, Yuri Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons. 2009.
  * [NURBS] Les Piegl & Wayne Tiller: The NURBS Book, Springer-Verlag 1995–1997 (2nd ed.).
  * [BE] Michael J. Borden, Michael A. Scott, John A. Evans, and Thomas J.R. Hughes: Isogeometric Finite Element Data Structures based on Bezier Extraction of NURBS, Int. J. Numer. Meth. Engng., 87: 15–47. doi: 10.1002/nme.2968, 2011.
  * [PetIGA] N. Collier, L. Dalcin, V.M. Calo: PetIGA: High-Performance Isogeometric Analysis, arxiv 1305.4452, 2013, http://arxiv.org/abs/1305.4452.
  * [SfePy] Robert Cimrman: SfePy - Write Your Own FE Application, arxiv 1404.6391, 2014, http://arxiv.org/abs/1404.6391.