Nathan Collier  committed faf283e

work on 4

  • Participants
  • Parent commits ba1f02d
  • Branches master

Comments (0)

Files changed (1)

File TUTORIAL4.rst

 .. role:: file(literal)
-Tutorial 4 - Disecting the Poission Codes
+Tutorial 4 - Dissecting the Poission Codes
 **Objective:** At the end of this tutorial you will understand 
-**Assumptions:** I assume 
+**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 most PetIGA applications, you only need to include
+:file:`petiga.h`. Including this file will also automatically include
+ #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. 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 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 right hand side.
+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. 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
+ #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;
+ }
+ #undef __FUNCT__
+ #define __FUNCT__ "main"
+ int main(int argc, char *argv[]) {
+   PetscErrorCode  ierr;
+   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+  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);
+  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);    
+  }
+  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);
+  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);
+  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
+  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;
+ }
 Higher Dimensions and Dimension Independent