include("../JuliaIGA.jl")
using Plots
plotly()
#pyplot()
Plots.PlotlyBackend()
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}$.
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.
JuliaIGA.BezierPlot(C,0.01)
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.
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.
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]
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.
ξ=[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)
ξ=[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
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)$$
#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)
We show how to plot a NURBS surface, in particular a torus. In particular we implemented two different algorithms to evaluate the NURBS surface.
We can tell NURBSPlot which option we prefer using the opt argument, by default we use the Projection algorithm.
ξ=[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
We first define the geometry where we are going to solve the Poisson problem.
#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)
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)
B=[
[0.0 0.0 0.0] [1.0 0.0 0.0];
[0.0 1.0 0.0] [1.0 1.0 0.0];
]
S = JuliaIGA.BezierSurface(1,1,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)
println(S.B)
B=[
[0.0 0.0 0.0] [0.5 0.0 0.0] [1.0 -0.5 0.0];
[0.0 1.0 0.0] [0.5 1.0 0.0] [1.0 0.5 0.0];
]
S = JuliaIGA.BezierSurface(1,2,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)