In [1]:

```
include("../JuliaIGA.jl")
using Plots
plotly()
#pyplot()
```

Out[1]:

Different function to draw curve and surface were implemented, in particular we implemented separate function to define and draw Bezier curves and surfaces, B-Spline curve and surfaces and last but not least we implemented as well a set of function to define and draw NURBS curves and surfaces.

We define the **Bernstein polynomial** of order $k$ and degree $n$ the following polynomial defined over $[0,1]$:
$$b^n_k(t)=\begin{pmatrix}n\\k\end{pmatrix}(1-t)^{n-k}t^k$$
we can then define the n-th **Bernstein polynomial** associated with a continuous function $f:[0,1]\to \mathbb{R}$ the polynomial:
$$p_n(t)=\sum_{k=0}^n f\Big(\frac{k}{n}\Big)b^n_k(t)$$
We define a Bezier curve $\vec{C}:[0,1]\to \mathbb{R}$, of control points $\{\vec{p}_0,\dots,\vec{p}_n\}$, the following function:
$$\vec{C}(t) = \sum_{k=0}^p \vec{p}_k b^n_k(t)$$
We used the De Casteljau algorithm to evaluate the Bezier curve in a given point.
In the following line of code we define a 2D Bezier curves of degree $2$ that has as control points: $\begin{pmatrix}0\\0\end{pmatrix}$,$\begin{pmatrix}0.5\\1\end{pmatrix}$,$\begin{pmatrix}1\\0\end{pmatrix}$.

In [2]:

```
C = JuliaIGA.BezierCurve2D(2,[[0 0],[0.5 1],[1 0]])
#We can access the information of the Bezier curve, after its definition.
println("Bezier curve of order",C.p,"and control points:")
println(C.V)
```

We also implemented the special function that use Plots to display the Bezier curves.

In [3]:

```
JuliaIGA.BezierPlot(C,0.01)
```

Out[3]:

We define instead a Bezier surface as the tensor product of two Bezier curves, i.e. $$ \vec{S}(s,t)=\sum_{j=0}^p\sum_{k=0}^q \vec{p}_{jk} b^q_j(s)b^p_k(t)$$ we see next how to define and plot a Bezier surface.

In [4]:

```
B=[
[0.0 0.0 0.0] [-5.0 18.0 5.0] [-7.0 32.0 7.0] [0.0 50.0 -3.0];
[12.0 0.0 -5.0] [13.0 17.0 -3.0] [14.0 31.0 10.0] [12.0 55.0 5.0];
[28.0 1.0 4.0] [29.0 15.0 3.0] [28.0 28.0 0.0] [29.0 45.0 -7.0];
[40.0 0.0 1.0] [43.0 18.0 -1.0] [37.0 32.0 5.0] [40.0 50.0 0.0];
]
S = JuliaIGA.BezierSurface(3,3,B)#We here define the following a bicubic Bezier surface of control points B.
println("Degere of the first Bezier curve, ",S.p,", degree of the second Bezier curve, ",S.q,".")
X,Y,Z = JuliaIGA.BezierPlot(S,0.01)#We here export the variable to build the graph of the surface.
surface(X,Y,Z)#We use the Plots library to plot the Bezier surface.
#JuliaIGA.ControlPlot(S)
```

Out[4]:

In [5]:

```
println(S.B)#Printing Control Points Matrix For The Bezier Surface
```

We can define the i-th B-spline base function of order $p$ on the domain $[0,1]$, associated with a knot vector $[\xi_0,\dots,\xi_{n+p+1}]$:
$$N_{i,0}(\xi) = \begin{cases}1\quad \xi_i \leq \xi \leq \xi_{i+1}\\ 0\end{cases}$$
$$N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p}-\xi_i}N_{i,p-1}(\xi) + \frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}N_{i+1,p-1}$$
we define a **piecewise polynomial B-Spline curve** of control points $\{\vec{B}_0,\dots,\vec{B}_n\}$ the following function $\vec{C}:[\xi_0,\xi_n]\to \mathbb{R}$:
$$C(\xi) = \sum_{i=1}^n \vec{B}_iN_{i,p}(\xi)$$
We used the deBoor recursive formula to evaluate the BSpline in a given set of points.

In [6]:

```
Î¾=[0.0,0.0,0.0,1.0,2.0,3.0,3.0,3.0]
V=[0.0,0.0,1.0,0.0,0.0]
C=JuliaIGA.BSplineCurve(2,V,Î¾)
JuliaIGA.BSplinePlot(C,0.01)
```

Out[6]:

In [7]:

```
Î¾=[0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0] # We define the knot vector
V=[[-4.0 -4.0],[-2.0 4.0],[2.0 -4.0],[4.0 4.0]] # We define the control points of the BSpline
C2 = JuliaIGA.BSplineCurve2D(3,V,Î¾) # We build the BSpline having Î¾ as knot vector and V as control polygon.
JuliaIGA.BSplinePlot(C2,0.01)# We plot the BSpline preaviusly defined
```

Out[7]:

We can describe using piece wise polynomial some shape such as a circle, to prevent this problem from occurring we introduce a particular mathematical object known as NURBS. A NURBS is the weighted-projection of a B-spline living in $\mathbb{R}^{d+1}$ on the $\mathbb{R}^d$ space. In particular given a knot vector $\xi$, a weight vector $\omega$ we define the associated NURBS base: $$R^p_i(\xi)=\frac{N_{i,p}(\xi)\omega_i}{\sum_{i=1}^n N_{i,p}(\xi)\omega_i}$$ where the $N_{i,p}(\xi)$ is the n-th B-spline base of order $p$. Given the NURBS base we are now able to define a NURBS curve of control points $\{\vec{B}_0,\dots,\vec{B}_n\}$: $$\vec{C}(\xi) = \sum{n}_{i=1} \vec{B}_i R^p_i(\xi)$$

In [8]:

```
#We define the knot vector
Î¾=[0.0,0.0,0.0,0.5*Ï€,0.5*Ï€,Ï€,Ï€,1.5*Ï€,1.5*Ï€,2*Ï€,2*Ï€,2*Ï€];
#We define the weight vector
Ï‰=[1.0,0.5*sqrt(2),1.0,0.5*sqrt(2),1.0,0.5*sqrt(2),1.0,0.5*sqrt(2),1.0];
#We define the control points
B=[[1.0 0.0],[1.0 1.0],[0.0 1.0],[-1.0 1.0],[-1.0 0.0],[-1.0 -1.0],[0.0 -1.0],[1.0 -1.0],[1.0 0.0]];
#We define a second order NURBS, of knot vector Î¾, weight vector Ï‰ and control points B
S = JuliaIGA.NURBSCurve2D(2,B,Ï‰,Î¾);
#We can now plot the NURBS Curve.
JuliaIGA.NURBSPlot(S,0.01)
```

Out[8]:

We show how to plot a NURBS surface, in particular a torus. In particular we implemented two different algorithms to evaluate the NURBS surface.

**Projection**we project the 4D Bspline surface associated with the NURBS in the 3D space.**Direct**we directly compute the NURBS from its representation using as a base the tensor product of $R_{i}^{p}(\xi)$. We can tell NURBSPlot which option we prefer using the opt argument, by default we use the Projection algorithm.

In [9]:

```
Î¾=[0,0,0,1,1,2,2,3,3,4,4,4];
Ï„=[0,0,0,1,1,2,2,3,3,4,4,4];
B=[
[5 0 -1] [6 0 -1] [6 0 0] [6 0 1] [5 0 1] [4 0 1] [4 0 0] [4 0 -1] [5 0 -1];
[5 5 -1] [6 6 -1] [6 6 0] [6 6 1] [5 5 1] [4 4 1] [4 4 0] [4 4 -1] [5 5 -1];
[0 5 -1] [0 6 -1] [0 6 0] [0 6 1] [0 5 1] [0 4 1] [0 4 0] [0 4 -1] [0 5 -1];
#
[-5 5 -1] [-6 6 -1] [-6 6 0] [-6 6 1] [-5 5 1] [-4 4 1] [-4 4 0] [-4 4 -1] [-5 5 -1];
[-5 0 -1] [-6 0 -1] [-6 0 0] [-6 0 1] [-5 0 1] [-4 0 1] [-4 0 0] [-4 0 -1] [-5 0 -1];
[-5 -5 -1] [-6 -6 -1] [-6 -6 0] [-6 -6 1] [-5 -5 1] [-4 -4 1] [-4 -4 0] [-4 -4 -1] [-5 -5 -1];
#
[0 -5 -1] [0 -6 -1] [0 -6 0] [0 -6 1] [0 -5 1] [0 -4 1] [0 -4 0] [0 -4 -1] [0 -5 -1];
[5 -5 -1] [6 -6 -1] [6 -6 0] [6 -6 1] [5 -5 1] [4 -4 1] [4 -4 0] [4 -4 -1] [5 -5 -1];
[5 0 -1] [6 0 -1] [6 0 0] [6 0 1] [5 0 1] [4 0 1] [4 0 0] [4 0 -1] [5 0 -1];
];
Ï‰=[
1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1;
1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2);
1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1;
#
1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2);
1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1;
1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2);
#
1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1;
1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2) 0.5 1/sqrt(2);
1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1 1/sqrt(2) 1;
];
S = JuliaIGA.NURBSurface(2,2,Î¾,Ï„,B,Ï‰);
println("NURBS Surface")
X,Y,Z = JuliaIGA.NURBSPlot(S,0.1,opt="Projection")
surface(X,Y,Z)
```

Out[9]:

We first define the geometry where we are going to solve the Poisson problem.

In [14]:

```
#include("../JuliaIGA.jl")
Î¾=[0,0,0,(1/3),(2/3),1,1,1];
Î·=[0,0,0,(1/3),(2/3),1,1,1];
B=[
[0. 1. 0] [0.2612 1.0 0] [0.7346 0.7346 0] [1.0 0.2612 0] [1.0 0. 0];
[0.0 1.25 0] [0.3265 1.25 0] [0.9182 0.9182 0] [1.25 0.3265 0] [1.25 0.0 0];
[0.0 1.75 0] [0.4571 1.75 0] [1.2856 1.2856 0] [1.75 0.4571 0] [1.75 0.0 0];
[0.0 2.25 0] [0.5877 2.25 0] [1.6528 1.6528 0] [2.25 0.5877 0] [2.25 0.0 0];
[0.0 2.5 0] [0.6530 2.5 0] [1.8365 1.8365 0] [2.5 0.6530 0] [2.5 0.0 0];
];
Ï‰=[
1.0 0.9024 0.8373 0.9024 1.0;
1.0 0.9024 0.8373 0.9024 1.0;
1.0 0.9024 0.8373 0.9024 1.0;
1.0 0.9024 0.8373 0.9024 1.0;
1.0 0.9024 0.8373 0.9024 1.0;
];
S = JuliaIGA.NURBSurface(2,2,Î¾,Î·,B,Ï‰);
X,Y,Z = JuliaIGA.NURBSPlot(S,0.1,opt="Projection")
surface(X,Y,Z)
```

Out[14]:

In [16]:

```
BS = JuliaIGA.BezierExtraction(S)
X,Y,Z = JuliaIGA.BezierPlot(BS,0.01)#We here export the variable to build the graph of the surface.
surface(X,Y,Z)#We use the Plots library to plot the Bezier surface.
JuliaIGA.ControlPlot(BS)
```

Out[16]: