-Creating Geometry with igakit (*)

-=================================

+Creating Geometry with igakit

+=============================

-**Objective:** At the end of this tutorial you will understand how

- geometries may be generated with igakit and then be read and used in PetIGA.

+**Objective:** At the end of this tutorial you will understand how geometries may be generated with igakit and then be read and used in PetIGA.

-**Assumptions:** I assume

+**Assumptions:** I assume that you are familiar with python and have both packages installed and working.

A Brief Overview of igakit

--------------------------

+We created igakit to enable both a interface into many low-level

+routines which are needed to manipulate knot vectors or control point

+grids as well as a high level interface to simplify the creation of

+NURBS geometries by hand. These geometries then can be read in to

+PetIGA for use as a mapping. We also use the package for

+post-processing solution results from PetIGA. Below we briefly sketch

+some of the objects available in the package.

+* :file:`igalib` - a python module which contains several

+ submodules. This is a python wrapper around the Fortran module

+ implemented in :file:`igalib.f90`.

+ * :file:`bsp` - contains interfaces for B-spline related functions such as

+ :file:`FindSpan`, :file:`EvalBasisFunsDers`, :file:`Greville`,

+ :file:`InsertKnot`, :file:`RemoveKnot`, :file:`DegreeElevate`

+ * :file:`iga` - contains interfaces for a few functions that can be

+ used to solve B-spline-based Galerkin problems in python

+ * :file:`KnotVector` - will create a uniform knot vector of

+ user-specified number of elements, polynomial order,

+ inter-element continuity

+ * :file:`BasisData` - returns a series of objects needing to

+* :file:`nurbs` - contains a dimension independent python class

+ (:file:`NURBS`) which simplifies the usage of routines contained in

+ the :file:`igalib` module to create and manipulate NURBS curves,

+* :file:`io` - a module which contains two python classes:

+ :file:`PetIGA` and :file:`VTK` which are used to create input/output

+ files suitable to read in by other software.

+* :file:`cad` - a python module which is a collection of python

+ functions which generate commonly used geometries (lines,circles)

+ and perform geometry operations (extrude, join, revolve, ruled,

+* :file:`plot` - a python module which enables plotting natively in

+In this tutorial we will not exhaust the possibilities of igakit for

+generating geometry, but rather explain how igakit geometries can be

+created and read into PetIGA. We encourage the reader to explore the

+demo directories of the igakit package for more samples in how igakit

+can be used to generate NURBS geometries.

+As an example, we will create the ever-popular patch test of

+isogeometric analysis: the quarter annulus. I will write this tutorial

+as if you are following me in the python interpreter. If you aren't

+using ipython_, I highly recommend it. First we import the circle

+functionality from the :file:`cad` module::

+ >>> from igakit.cad import circle

+Most functions we have implement come with documentation (of various

+degrees of completeness) built-in to the source code. You can view

+where you will see something like the following::

+ circle(radius=1, center=None, angle=None)

+ Construct a NURBS circular arc or full circle

+ radius : float, optional

+ center : array_like, optional

+ angle : float or 2-tuple of floats, optional

+The third argument of the circle is of interest here. The angle lets

+us choose a sweep to define our annulus. We can use this to create a

+ >>> from numpy import pi

+ >>> c1 = circle(angle=(0,pi/2.))

+The circle function returns a NURBS object. To get an idea of what you

+can possibly do with this object, examine the help for it::

+You will see a long list of functions and examples of attributes and

+operations which are possible in the NURBS class. For now, let's

+inspect our curve to see what was created. If you are using ipython_,

+and then hit the TAB key. You will see a more concise list of what can

+ c1.array c1.control c1.evaluate c1.insert c1.refine c1.scale c1.transform

+ c1.boundary c1.copy c1.extract c1.knots c1.remap c1.shape c1.translate

+ c1.breaks c1.degree c1.fields c1.move c1.remove c1.slice c1.transpose

+ c1.clamp c1.dim c1.gradient c1.plot c1.reverse c1.spans c1.unclamp

+ c1.clone c1.elevate c1.greville c1.points c1.rotate c1.swap c1.weights

+So we see that we can inspect the knot vector::

+ (array([ 0., 0., 0., 1., 1., 1.]),)

+As expected, we created a single span of quadratics. The control point

+grid can also be inspected::

+ array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],

+ [ 7.07106781e-01, 7.07106781e-01, 0.00000000e+00, 7.07106781e-01],

+ [ 2.22044605e-16, 1.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

+which are the control points in the 4D space. The first three columns

+are the x-y-z locations of the control points and the last are the

+weights. Now we create a second circular arc, this time with a radius

+ >>> c2 = circle(radius=2,angle=(0,pi/2.))

+We also have some plotting functionality built into igakit. We will

+use this now to visually check that our curves have been created as

+ >>> from igakit.plot import plt

+ >>> plt.plot(c1,color='b')

+ >>> plt.plot(c2,color='g')

+You will see a 3D window open with the two curves we created that

+should look something like the following image.

+.. image:: ./figs/two_curves.png

+Now, we will create a `ruled surface`_ using another function in the

+:file:`cad` module, and then also check the result with a plot::

+ >>> from igakit.cad import ruled

+.. image:: ./figs/qannulus.png

Reading Geometry Into PetIGA

----------------------------

+Now while we have the geometry in hand, we need to prepare it for

+analysis. Let's inspect the polynomial orders of the surface that we

+Here we see that the first direction is the quadratic curve that we

+had originally, but the ruled operation linearly blended the bounding

+curves together. Let's assume that we prefer to use equal polynomial

+orders in each direction for the annulus. So we need to degree elevate

+the second direction by one.

+We can also inspect the knot vectors::

+ (array([ 0., 0., 0., 1., 1., 1.]), array([ 0., 0., 0., 1., 1., 1.]))

+Where we see that we have a single nonzero span of

+bi-quadratics. Since PetIGA does not handle geometry modification

+internally we need to insert knots here such that the function space

+matches the one we will use in the analysis portion. To this end, we

+will create a 10 by 10 element mesh of *C1* quadratics by knot

+ >>> from numpy import linspace

+ >>> to_insert = linspace(0,1,11)[1:-1]

+ >>> srf.refine(0,to_insert)

+ >>> srf.refine(1,to_insert)

+ (array([ 0. , 0. , 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1. , 1. ]),

+ array([ 0. , 0. , 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1. , 1. ]))

+Now the geometry is ready for analysis and to be exported to

+PetIGA. We use functionality from the :file:`io` module of igakit.

+ >>> from igakit.io import PetIGA

+ >>> PetIGA().write("quarter_annulus.dat",srf)

+This will create a datafile that can be read by PetIGA. Now we will

+use this geometry in the Poisson code with the setup from the previous

+tutorial, :ref:`TUTORIAL2`. Exit python and compile the demo code

+:file:`$PETIGA_DIR/demo/Poisson.c`. I am assuming that the file

+:file:`quarter_annulus.dat` that we just created is in this directory. If it

+is not here, then you will need to move it here. Then run the PetIGA

+code with these options::

+ Poisson -iga_geometry quarter_annulus.dat -iga_view -save

+This code solves the Poisson equation with unit Dirichlet boundary

+conditions everywhere and a unit body force. When we execute this

+code, we can see what function space was used. Note that this

+corresponds to what we created in igakit::

+ IGA: dim=2 dof=1 order=2 geometry=2 rational=1 property=0

+ Axis 0: periodic=0 degree=2 quadrature=3 processors=1 nodes=12 elements=10

+ Axis 1: periodic=0 degree=2 quadrature=3 processors=1 nodes=12 elements=10

+ Partitioning - nodes: sum=144 min=144 max=144 max/min=1

+ Partitioning - elements: sum=100 min=100 max=100 max/min=1

+Unlike the previous tutorial, the geometry is now flagged as 2D and

+the rational is flagged to True. If the input geometry is a NURBS then

+PetIGA detects and compensates for this automatically. Also, if a

+geometry is used, the basis is mapped automatically with zero code

+changes required by the user.

+Just as in the previous tutorial, the :file:`-save` option generates

+two files (:file:`PoissonGeometry.dat` and

+:file:`PoissonSolution.dat`). Then we can run the same visualization

+which generates :file:`PoissonVTK.vtk`. However, this time the

+resulting VTK file encodes the solution on the mapped geometry as in

+.. image:: ./figs/quarter_annulus_poisson.png

+Now, what if we did not want to take advantage of the quarter

+symmetry? How do we generate this geometry? Here I will write a single

+script which extends the above ideas to all four quadrants::

+ from igakit.io import PetIGA

+ from numpy import linspace

+ c2 = cad.circle(radius=2)

+ srf = cad.ruled(c1,c2).elevate(1,1)

+ to_insert = linspace(0,0.25,11)[1:-1]

+ srf.refine(0,i*0.25+to_insert)

+ srf.refine(1,i*0.25+to_insert)

+ PetIGA().write("annulus.dat",srf)

+But if run this geometry through the Poisson solver, this geometry

+will not produce the desired effect, but rather as shown in the

+.. image:: ./figs/annulus_no.png

+This is because we set all boundary conditions to one. We need to set

+the boundaries in the first parametric direction to have periodic

+boundary conditions. We need to add a line to our script, just before

+we output the geometry using the PetIGA class::

+ from igakit.io import PetIGA

+ from numpy import linspace

+ c2 = cad.circle(radius=2)

+ srf = cad.ruled(c1,c2).elevate(1,1)

+ to_insert = linspace(0,0.25,11)[1:-1]

+ srf.refine(0,i*0.25+to_insert)

+ srf.refine(1,i*0.25+to_insert)

+ srf.unclamp(0) # <----------------------- added this line

+ PetIGA().write("annulus.dat",srf)

+This will modify the geometry to use unclamped knot vectors as

+opposed to clamped (or open) ones. Now we can rerun the code,

+specifying the periodicity::

+ ./Poisson -iga_geometry annulus.dat -iga_view -save -iga_periodic 1,0

+ IGA: dim=2 dof=1 order=2 geometry=2 rational=1 property=0

+ Axis 0: periodic=1 degree=2 quadrature=3 processors=1 nodes=44 elements=40

+ Axis 1: periodic=0 degree=2 quadrature=3 processors=1 nodes=39 elements=37

+ Partitioning - nodes: sum=1716 min=1716 max=1716 max/min=1

+ Partitioning - elements: sum=1480 min=1480 max=1480 max/min=1

+You will see now that the first direction is flagged periodic and the

+resulting plot now displays the proper symmetry.

+.. image:: ./figs/annulus_yes.png

+.. _ipython: http://ipython.org/

+.. _ruled surface: http://en.wikipedia.org/wiki/Ruled_surface