# 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.

**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.

`igalib`

- a python module which contains several submodules. This is a python wrapper around the Fortran module implemented in`igalib.f90`

.`bsp`

- contains interfaces for B-spline related functions such as`FindSpan`

,`EvalBasisFunsDers`

,`Greville`

,`InsertKnot`

,`RemoveKnot`

,`DegreeElevate`

`iga`

- contains interfaces for a few functions that can be used to solve B-spline-based Galerkin problems in python`KnotVector`

- will create a uniform knot vector of user-specified number of elements, polynomial order, inter-element continuity`BasisData`

- returns a series of objects needing to integrate weak form

`nurbs`

- contains a dimension independent python class (`NURBS`

) which simplifies the usage of routines contained in the`igalib`

module to create and manipulate NURBS curves, surfaces, and volumes.`io`

- a module which contains two python classes:`PetIGA`

and`VTK`

which are used to create input/output files suitable to read in by other software.`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, sweep).`plot`

- a python module which enables plotting natively in python.

## Creating a Geometry¶

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 `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 this by:

```
>>> help(circle)
```

where you will see something like the following:

```
circle(radius=1, center=None, angle=None)
Construct a NURBS circular arc or full circle
Parameters
----------
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 circular arc:

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

```
>>> help(c1)
```

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, you can type:

```
>>> c1.
```

and then hit the TAB key. You will see a more concise list of what can be done to this curve:

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

```
>>> c1.knots
(array([ 0., 0., 0., 1., 1., 1.]),)
```

As expected, we created a single span of quadratics. The control point grid can also be inspected:

```
>>> c1.control
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 of 2:

```
>>> 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 expected:

```
>>> from igakit.plot import plt
>>> plt.plot(c1,color='b')
>>> plt.plot(c2,color='g')
>>> plt.show()
```

You will see a 3D window open with the two curves we created that should look something like the following image.

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

module, and then also check the result with a plot:

```
>>> from igakit.cad import ruled
>>> srf = ruled(c1,c2).transpose()
>>> plt.plot(srf)
>>> plt.show()
```

## 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 created:

```
>>> srf.degree
(1,2)
```

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.

```
>>> srf.elevate(0,1)
>>> srf.degree
(2,2)
```

We can also inspect the knot vectors:

```
>>> srf.knots
(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
insertion:

```
>>> from numpy import linspace
>>> to_insert = linspace(0,1,11)[1:-1]
>>> srf.refine(0,to_insert)
>>> srf.refine(1,to_insert)
>>> srf.knots
(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 `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, Post-processing Results with igakit. Exit python and compile the demo code
`$PETIGA_DIR/demo/Poisson.c`

. I am assuming that the 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 `-save`

option generates
two files (`PoissonGeometry.dat`

and
`PoissonSolution.dat`

). Then we can run the same visualization
script:

```
python PoissonPlot.py
```

which generates `PoissonVTK.vtk`

. However, this time the
resulting VTK file encodes the solution on the mapped geometry as in
the following image.

## Periodicity¶

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 import cad
from igakit.io import PetIGA
from numpy import linspace
c1 = cad.circle()
c2 = cad.circle(radius=2)
srf = cad.ruled(c1,c2).transpose().elevate(0,1)
to_insert = linspace(0,0.25,11)[1:-1]
for i in range(4):
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 following figure.

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 import cad
from igakit.io import PetIGA
from numpy import linspace
c1 = cad.circle()
c2 = cad.circle(radius=2)
srf = cad.ruled(c1,c2).transpose().elevate(0,1)
to_insert = linspace(0,0.25,11)[1:-1]
for i in range(4):
srf.refine(0,i*0.25+to_insert)
srf.refine(1,i*0.25+to_insert)
srf.unclamp(1) # <----------------------- 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.