.. role:: option(literal)
.. role:: file(literal)
.. _TUTORIAL4:

Dissecting the Poisson Codes
=============================

**Objective:** At the end of this tutorial you will understand how a simple linear solver works in PetIGA.

**Assumptions:** I assume that you are familiar with C and in particular understand pointers_. 

Detailed Look at :file:`Poisson1D.c`
------------------------------------

In this tutorial we will look at the Poisson codes and explain each
portion in more detail. It is my hope that in explaining how the demo
codes work, this will help you see how you can use PetIGA to address
your own problems. In this demo I will be discussing this_ demo code
in sections. It might be useful to first scan through this code to
familiarize yourself with its organization.

In most PetIGA applications, you only need to include
:file:`petiga.h`. Including this file will also automatically include
:file:`petsc.h`.

::

 #include "petiga.h"

The first function we encounter is the :file:`System` function. This
function is the key to understanding how to use PetIGA to solve your
own problems. For a linear, steady problem, this function represents
the evaluation of the bilinear and linear form at a quadrature point.

At a glance, a few things may strike you as strange particularly if
you are not familiar with PETSc. First, you will note the
:file:`#undef` and :file:`#define` lines before the function
declaration. This is a mechanism that PETSc uses such that when its
functions fail, the entire call stack is returned to the user. This is
typically something that you only get when running the
debugger. Having it always return on error reduces the time it takes
you to locate errors.

You will also see variable types :file:`PetscReal` and
:file:`PetscScalar` instead of the expected :file:`double`. These are
PETSc types which for most applications are actually of the type
:file:`double`. However, you can configure PETSc to use, for example,
quadruple precision or complex numbers. In this case, using PETSc
types means that your code automatically ports to these situations
with no additional work for you.

The function has as its arguments a :file:`IGAPoint`, two PetscScalars
and a void pointer. For now, ignore the existence of the void
pointer. The :file:`IGAPoint p` is the quadrature/collocation
point. We have built in all the information that you will need to
evaluate the bilinear and linear form at this point. The two
PetscScalar pointers are for the discretized bilnear and linear form
respectively. The variable :file:`K` can be thought of as a flattened
point contribution to the element stiffness matrix. The variable
:file:`F` is the point contribution to the element load vector.

The basis functions are part of the IGAPoint :file:`p` and so we need
to simply assign some pointers to their location. The pointers
:file:`N0` and :file:`N1` represent the 0th and 1st order derivative
of the basis functions respectively. The pointers used and their
casting is a bit technical and ugly, a consequence of using the C
language. However, we emphasize that it is sufficient to copy these
lines into your own code and understand what it is the pointers
address. If a geometry is used, then the basis is mapped and/or made
rational if using NURBS. Otherwise, these are just the B-spline basis
functions in the parametric space. Your code does not need to know the
difference, it is all handled internally.

Then we simply code a double loop over the local number of basis
functions, :file:`nen`. As a convention, we use :file:`a` to loop over
the test functions and :file:`b` to loop over the trial functions. A
dot product of the basis function derivatives is assembled into the
stiffness matrix and the unit body force is assembled into the load
vector. The quadrature weights and potential Jacobian of the mapping
are applied internally.

::

 #undef  __FUNCT__
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   const PetscReal *N0,(*N1)[1];
   IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
   IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
     PetscReal Na_x = N1[a][0];
     for (b=0; b<nen; b++) {
       PetscReal Nb_x = N1[b][0];
       K[a*nen+b] = Na_x * Nb_x;
     }
     F[a] = Na * 1.0;
   }
   return 0;
 }

Below this we code the program's main function. As in all PETSc codes,
we intially call an initialization routine. Consider this a
rite-of-passage which is required for PETSc's internals and
parallelism. However, another oddity may jump out at you. You will see
that we have defined an error code :file:`ierr` and then after the
PETSc call, we postpend a function call :file:`CHKERRQ(ierr)`. If you
glance at the rest of the code, you will see this pattern appears at
most every line of the code (yes they are ugly, be patient, eventually
you won't notice them). Every function in PETSc/PetIGA returns such a
code which is really just an integer whose value corresponds to a
error type. The function call at the end is really a macro which will
cause the application to fail and report the proper error
messages/stack.

::

 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {

   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);

Now we create an :file:`IGA`, your one-stop-shop for the
IGA-discretization. We set the dimensions of the space as well as the
number of degrees of freedom per basis. Basic usuage of an :file:`IGA`
in PetIGA assume that the same function space is used for all
test/trial spaces in the problem. This is a constraint that we will
later show how to relax. The :file:`XXXSetFromOptions` calls in PETSc
are very important. This call is what populates the list of options
that you see when you run the code with the :file:`-help` option. The
final call is an internal setup of many data structures which we need
to compute. This is considered a final step of creating the
:file:`IGA`. There are many other ways to initialize an :file:`IGA`,
this way is just the simplest and uses the default discretizations.

::

  IGA iga;
  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
  ierr = IGASetDim(iga,1);CHKERRQ(ierr);
  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
  ierr = IGASetUp(iga);CHKERRQ(ierr);

In this block we set the boundary conditions. Again, there are many
ways you can set boundary conditions in PetIGA. Here we will set a
single value on all basis functions on the boundaries. The
:file:`IGABoundary` represents a side in 1D, a boundary edge in 2D,
and a boundary face in 3D. You get the boundary object from the
:file:`IGA` itself by specifying a direction (which parametric
dimension) and then a side (0 or 1 reflecting the min and max of the
axis). With the :file:`IGABoundary` in hand, then you can set a
constant value on all associated basis functions of a particular
component of the solution. Since this is a scalar problem, we set the
0th component.

::

  IGABoundary bnd;
  PetscInt dir=0,side;
  PetscScalar value = 1.0;
  for (side=0; side<2; side++) {
    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
    ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);    
  }

Now we setup and compute the linear system of our problem,
:file:`Ax=b`. Both :file:`Mat` and :file:`Vec` are PETSc types which
encapsulate sparsity and parallel communication. You do not need to
understand how these work internally--they are given to you by the
:file:`IGA`. The call to :file:`IGASetUserSystem` associates the
:file:`System` function we wrote earlier to the discretization. The
:file:`IGA` uses this function in the next call to assemble the linear
system.

::

  Mat A;
  Vec x,b;
  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);

Then we need a :file:`KSP` to solve the linear system. This we also
get from our :file:`IGA`. This merely ensures that the :file:`KSP` and
the :file:`IGA` use the same communicator (a detail of parallelism
which you can ignore). Then we need to set the problem operators by
KSPSetOperators_. The first assignment of :file:`A` is for the problem
operator. The second assignment of :file:`A` specifies the matrix from
which PETSc will create a preconditioner. In this case, and in many
cases, we want the preconditioner to be created from original system
matrix. We also specify a flag for the :file:`KSP` that our
preconditioner has the same nonzero pattern as the matrix
:file:`A`. As before with the :file:`IGA`, we add the
:file:`XXXSetFromOptions` call which will allows us to control
everything about the :file:`KSP` from the commandline. Finally, we
solve.

::
  
  KSP ksp;
  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

This call to file:`VecView` is optional. It is what will draw a rough
plot of the degrees of freedom, as described in :ref:`TUTORIAL2`

::

  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);

Finally, we destroy the objects that we created and issue a finalize
call to PETSc.

::
  
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = IGADestroy(&iga);CHKERRQ(ierr);
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
 }

That is all that you need to do, everything else in handled internally
in the PetIGA library. We feel that our approach is beneficial because
it leads to shorter codes which focus on the physics of the
problem. We also emphasize that this code will work in serial or in
parallel without the author needing to know anything about parallelism
aside from having to initialize objects with communicators.

Higher Dimensions and Dimension Independent
-------------------------------------------

How different is the code when changing dimensions of the problem?
There is very little which depends on the problem dimension in the
application code. If we compare to
:file:`PETIGA_DIR/demo/Poisson2d.c`, we have a few lines in our
:file:`System` function which are different. The pointers to the shape
function gradients are now 2D::

   const PetscReal *N0,(*N1)[1];

::

   const PetscReal *N0,(*N1)[2];

and the dot product is now 2D::

       K[a*nen+b] = Na_x * Nb_x;

::

       K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;

In :file:`main` we also have a few lines::

   ierr = IGASetDim(iga,1);CHKERRQ(ierr);

::

   ierr = IGASetDim(iga,2);CHKERRQ(ierr);

and we need to set more boundary conditions. In fact, this reflects
that a dimension independent code_ is simple to obtain. 

.. _pointers: http://en.wikipedia.org/wiki/Pointer_(computer_programming)
.. _this: https://bitbucket.org/dalcinl/petiga/src/f07e8f23a997aadf63f3809b1c5ac16aadb38307/demo/Poisson1D.c?at=default
.. _KSPSetOperators: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetOperators.html
.. _Tutorial 2: 
.. _code: https://bitbucket.org/dalcinl/petiga/src/f07e8f23a997aadf63f3809b1c5ac16aadb38307/demo/Poisson.c?at=default