JuliaIGA

In [1]:
include("../JuliaIGA.jl")
using Plots
plotly()
#pyplot()
Out[1]:
Plots.PlotlyBackend()

CAD

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.

Bezier Curves And Surface

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)
Bezier curve of order2and control points:
Array{Real,2}[[0.0 0.0], [0.5 1.0], [1.0 0.0]]

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)
Degere of the first Bezier curve, 3, degree of the second Bezier curve, 3.
Out[4]:
In [5]:
println(S.B)#Printing Control Points Matrix For The Bezier Surface
Real[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]

B-Spline Curves And 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]:

NURBS Curves And Surface

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.

  1. Projection we project the 4D Bspline surface associated with the NURBS in the 3D space.
  2. 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)
NURBS Surface
Out[9]:

Bezier Extraction

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]: