Commits

Lisandro Dalcin  committed 9129d59

More work in collocation and update tests

  • Participants
  • Parent commits 3ba4e3d

Comments (0)

Files changed (6)

File include/petiga.h

   PETSCHEADER(struct _IGAOps);
   PetscBool  setup;
   PetscInt   setupstage;
+  PetscBool  collocation;
 
   IGAUserOps userops;
   VecType    vectype;
   IGAAxis     axis[3];
   IGARule     rule[3];
   IGABoundary boundary[3][2];
+  IGABasis    basis[3];
 
-  IGABasis    basis[3];
   IGAElement  iterator;
 
   PetscInt    geometry;
   Vec        vwork[16];
   Vec        natural;
   VecScatter n2g,g2n;
-
-  /* stuff added for collocation */
-  PetscBool   collocation;
-  IGABasis    node_basis[3];
-  IGAElement  node_iterator;
 };
 
 PETSC_EXTERN PetscClassId IGA_CLASSID;

File src/petiga.c

   }
   ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
 
-  for (i=0; i<3; i++) /* collocation */
-    {ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);}
-  ierr = IGAElementDestroy(&iga->node_iterator);CHKERRQ(ierr);
-
   ierr = IGAReset(iga);CHKERRQ(ierr);
   ierr = PetscHeaderDestroy(&iga);CHKERRQ(ierr);
   PetscFunctionReturn(0);
   ierr = DMDestroy(&iga->node_dm);CHKERRQ(ierr);
 
   ierr = IGAElementReset(iga->iterator);CHKERRQ(ierr);
-  ierr = IGAElementReset(iga->node_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
       ierr = PetscViewerASCIIPrintf(viewer,"Partitioning - elements: sum=%D  min=%D  max=%D  max/min=%g\n",
                                     isum[1],imin[1],imax[1],(double)imax[1]/(double)imin[1]);CHKERRQ(ierr);
     }
-    if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
+    if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)
+      {
         PetscMPIInt rank; PetscInt *ranks = iga->proc_ranks;
         PetscInt *nnp = iga->node_lwidth, tnnp = 1, *snp = iga->node_lstart;
         PetscInt *nel = iga->elem_width,  tnel = 1, *sel = iga->elem_start;
         ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
         ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
       }
-    }
+  }
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "IGASetUseCollocation"
 PetscErrorCode IGASetUseCollocation(IGA iga,PetscBool collocation)
 {
-  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
     SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
             "Cannot change collocation after IGASetUp()");
   iga->collocation = collocation;
-  if (collocation && iga->setup) {
-    PetscInt i, dim = iga->dim;
-    PetscBool periodic = PETSC_FALSE;
-    for (i=0; i<dim; i++) if(iga->axis[i]->periodic) periodic = PETSC_TRUE;
-    if (periodic) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
-                          "Collocation not supported with periodicity");
-  }
-  if (collocation && iga->setup) {
-    PetscInt i;
-    for (i=0; i<3; i++) {
-      ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);
-      ierr = IGABasisCreate(&iga->node_basis[i]);CHKERRQ(ierr);
-      ierr = IGABasisInitCollocation(iga->node_basis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
-    }
-    ierr = IGAElementDestroy(&iga->node_iterator);CHKERRQ(ierr);
-    ierr = IGAElementCreate(&iga->node_iterator);CHKERRQ(ierr);
-    iga->node_iterator->collocation = PETSC_TRUE;
-    ierr = IGAElementInit(iga->node_iterator,iga);CHKERRQ(ierr);
-  }
   PetscFunctionReturn(0);
 }
 
       }
     }
 
+    /* Order */
+    if (order < 0) for (i=0; i<dim; i++) order = PetscMax(order,iga->axis[i]->p);
+    order = PetscMax(order,1); order = PetscMin(order,3);
+    ierr = PetscOptionsInt("-iga_order","Order","IGASetOrder",order,&order,&flg);CHKERRQ(ierr);
+    if (flg) { ierr = IGASetOrder(iga,order);CHKERRQ(ierr);}
+
     /* Quadrature */
     for (i=0; i<dim; i++) if (quadr[i] < 1) quadr[i] = iga->axis[i]->p + 1;
     ierr = PetscOptionsIntArray("-iga_quadrature","Quadrature points","IGARuleInit",quadr,(nq=dim,&nq),&flg);CHKERRQ(ierr);
         if (q > 0) {ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);}
       }
 
-    /* Order */
-    if (order < 0) for (i=0; i<dim; i++) order = PetscMax(order,iga->axis[i]->p);
-    order = PetscMax(order,1); order = PetscMin(order,3);
-    ierr = PetscOptionsInt("-iga_order","Order","IGASetOrder",order,&order,&flg);CHKERRQ(ierr);
-    if (flg) { ierr = IGASetOrder(iga,order);CHKERRQ(ierr);}
-
   setupcalled:
 
     /* Matrix and Vector type */
   iga->order = PetscMax(iga->order,1); /* XXX */
   iga->order = PetscMin(iga->order,3); /* XXX */
 
-  for (i=0; i<3; i++) {
-    if (i >= iga->dim) {ierr = IGARuleReset(iga->rule[i]);CHKERRQ(ierr);}
-    if (iga->rule[i]->nqp < 1) {ierr = IGARuleInit(iga->rule[i],iga->axis[i]->p + 1);CHKERRQ(ierr);}
-    ierr = IGABasisInitQuadrature(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
+  if (iga->collocation) {
+    for (i=0; i<iga->dim; i++)
+      if(iga->axis[i]->periodic)
+        SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,"Collocation not supported with periodicity");
   }
+
+  for (i=0; i<3; i++)
+    if (!iga->collocation && (i >= iga->dim || iga->rule[i]->nqp < 1))
+      {ierr = IGARuleInit(iga->rule[i],iga->axis[i]->p + 1);CHKERRQ(ierr);}
+    else
+      {ierr = IGARuleReset(iga->rule[i]);CHKERRQ(ierr);}
+
+  for (i=0; i<3; i++)
+    if (!iga->collocation)
+      {ierr = IGABasisInitQuadrature(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);}
+    else
+      {ierr = IGABasisInitCollocation(iga->basis[i],iga->axis[i],iga->order);CHKERRQ(ierr);}
+
   ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
 
   ierr = IGASetUp_View(iga);CHKERRQ(ierr);
-
-  if (iga->collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}  /* collocation */
-
   PetscFunctionReturn(0);
 }
 

File src/petigaelem.c

   IGACheckSetUp(iga,1);
   ierr = IGAElementReset(element);CHKERRQ(ierr);
   element->parent = iga;
+  element->collocation = iga->collocation;
 
   element->dof = iga->dof;
   element->dim = iga->dim;
   start = iga->elem_start;
   width = iga->elem_width;
   sizes = iga->elem_sizes;
-  BD = element->collocation ? iga->node_basis : iga->basis;
+  BD    = iga->basis;
 
   { /* */
     PetscInt i,dim = element->dim;
     PetscInt nel=1,nen=1,nqp=1;
-    for (i=0; i<3; i++)
+    for (i=0; i<3; i++) {
+      element->ID[i] = 0;
       element->BD[i] = BD[i];
+    }
     for (i=0; i<dim; i++) {
       element->start[i] = start[i];
       element->width[i] = width[i];
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidPointer(element,2);
   *element = iga->iterator;
-  if (PetscUnlikely(iga->collocation)) /* collocation */
-    *element = iga->node_iterator;
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "IGABeginElement"
 PetscErrorCode IGABeginElement(IGA iga,IGAElement *element)
 {
-  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidPointer(element,2);
   IGACheckSetUp(iga,1);
-  ierr = IGAGetElement(iga,element);CHKERRQ(ierr);
+  *element = iga->iterator;
   (*element)->index = -1;
   PetscFunctionReturn(0);
 }
 
 /* -------------------------------------------------------------------------- */
 
-static void AddFixa(IGAElement element,IGABoundary b,PetscInt a)
-{
-  if (b->count) {
-    PetscInt dof = element->dof;
-    PetscInt count = element->nfix;
-    PetscInt *index = element->ifix;
-    PetscScalar *value = element->vfix;
-    PetscInt j,k,n = b->count;
-    for (k=0; k<n; k++) {
-      PetscInt idx = a*dof + b->field[k];
-      PetscScalar val = b->value[k];
-      for (j=0; j<count; j++)
-        if (index[j] == idx) break;
-      if (j == count) count++;
-      index[j] = idx;
-      value[j] = val;
-    }
-    element->nfix = count;
-  }
-}
-
 EXTERN_C_BEGIN
 extern void IGA_BoundaryArea_2D(const PetscInt[],PetscInt,PetscInt,
                                 PetscInt,const PetscReal[],
   return A;
 }
 
+static void AddFixa(IGAElement element,IGABoundary b,PetscInt a)
+{
+  if (b->count) {
+    PetscInt dof = element->dof;
+    PetscInt count = element->nfix;
+    PetscInt *index = element->ifix;
+    PetscScalar *value = element->vfix;
+    PetscInt j,k,n = b->count;
+    for (k=0; k<n; k++) {
+      PetscInt idx = a*dof + b->field[k];
+      PetscScalar val = b->value[k];
+      for (j=0; j<count; j++)
+        if (index[j] == idx) break;
+      if (j == count) count++;
+      index[j] = idx;
+      value[j] = val;
+    }
+    element->nfix = count;
+  }
+}
+
 static void AddFlux(IGAElement element,IGABoundary b,PetscInt a,PetscReal A)
 {
   if (b->nload) {
       for (side=0; side<2; side++) {
         IGABoundary b = AtBoundary(element,dir,side);
         if(b && b->nload) {
-          PetscInt  f,n = b->nload;
-          PetscReal normal[3] = {0.0,0.0,0.0};
-          PetscReal *dshape;
+          PetscInt  f, n = b->nload;
+          PetscReal *dshape, normal[3] = {0.0,0.0,0.0};
           if (!element->geometry) {
             normal[dir] = side ? 1.0 : -1.0;
             dshape = element->basis[1];
             dshape = element->shape[1];
           }
           for (f=0; f<n; f++) {
-            PetscInt  c = b->iload[f];
-            PetscReal v = b->vload[f];
-            PetscInt  a,j;
+            PetscInt    c = b->iload[f];
+            PetscScalar v = b->vload[f];
+            PetscInt    a,j;
             for (j=0; j<N; j++) K[c*N+j] = 0.0;
             for (a=0; a<nen; a++)
               K[c*N+a*dof+c] = DOT(dim,&dshape[a*dim],normal);
           }
         }
         if (b && b->count) {
-          PetscInt  f,n = b->count;
-          PetscReal *shape;
-          if (!element->geometry) {
-            shape = element->basis[0];
-          } else {
-            shape = element->shape[0];
-          }
+          PetscInt  f, n = b->count;
+          PetscReal *shape = element->basis[0];
           for (f=0; f<n; f++) {
-            PetscInt  c = b->field[f];
-            PetscReal v = b->value[f];
-            PetscInt  a,j;
+            PetscInt    c = b->field[f];
+            PetscScalar v = b->value[f];
+            PetscInt    a,j;
             for (j=0; j<N; j++) K[c*N+j] = 0.0;
             for (a=0; a<nen; a++)
               K[c*N+a*dof+c] = shape[a];
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
-      if (element->rowmap[0] == element->colmap[a]) {
-        PetscInt    c = k % element->dof;
-        PetscScalar v = element->vfix[f];
-        PetscScalar u = element->ufix[f];
-        F[c] = u - v;
-      }
+      PetscInt c = k % element->dof;
+      if (element->rowmap[0] == element->colmap[a])
+        {
+          PetscScalar v = element->vfix[f];
+          PetscScalar u = element->ufix[f];
+          F[c] = u - v;
+        }
     }
   }
   PetscFunctionReturn(0);
     PetscInt f,n;
     n = element->nfix;
     for (f=0; f<n; f++) {
-      PetscInt k = element->ifix[f];
-      PetscInt i,j;
+      PetscInt i,j,k=element->ifix[f];
       for (i=0; i<M; i++) J[i*N+k] = 0.0;
       for (j=0; j<N; j++) J[k*N+j] = 0.0;
       J[k*N+k] = 1.0;
   PetscFunctionReturn(0);
  collocation:
   {
-    PetscInt M = element->neq * element->dof;
-    PetscInt N = element->nen * element->dof;
+    PetscInt nen = element->nen;
+    PetscInt dof = element->dof;
     PetscInt f,n;
     n = element->nfix;
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
-      PetscInt a = k / element->dof;
-      if (element->rowmap[0] == element->colmap[a]) {
-        PetscInt c = k % element->dof;
-        PetscInt i,j;
-        for (i=0; i<M; i++) J[i*N+k] = 0.0;
-        for (j=0; j<N; j++) J[c*N+j] = 0.0;
-        J[c*N+k] = 1.0;
-      }
+      PetscInt a = k / dof;
+      PetscInt c = k % dof;
+      if (element->rowmap[0] == element->colmap[a])
+        {
+          PetscInt  i,j,N=nen*dof;
+          PetscReal *shape = element->basis[0];
+          for (j=0; j<N; j++) J[c*N+j] = 0.0;
+          for (i=0; i<nen; i++)
+            J[c*N+i*dof+c] = shape[i];
+        }
     }
   }
   PetscFunctionReturn(0);

File src/petigamat.c

   PetscReal *U = iga->axis[dir]->U;
   PetscInt k;
 
+  if (PetscUnlikely(iga->collocation)) {
+    k = iga->basis[dir]->offset[i];
+    *first = k; *last  = k + p; return;
+  }
+
   /* compute index of the leftmost overlapping basis */
   k = i;
   while (U[k]==U[k+1]) k++; /* XXX Using "==" with floating point values ! */
     *first = k - s - n;
   }
 
-  if (PetscUnlikely(iga->collocation)) { /* collocation */
-    PetscInt offset = iga->node_basis[dir]->offset[i];
-    *first = offset;
-    *last  = offset + p;
-  }
 }
 
 PETSC_STATIC_INLINE

File test/Test_SNES_2D.c

   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "FunctionCollocation"
+PetscErrorCode FunctionCollocation(IGAPoint p,const PetscScalar *Ue,PetscScalar *Fe,void *ctx)
+{
+  PetscReal xy[2];
+  IGAPointFormPoint(p,xy);
+  PetscReal x = xy[0];
+  PetscReal y = xy[1];
+
+  PetscScalar u0[4],u1[4][2],u2[4][2][2];
+  IGAPointFormValue(p,Ue,&u0[0]);
+  IGAPointFormGrad (p,Ue,&u1[0][0]);
+  IGAPointFormHess (p,Ue,&u2[0][0][0]);
+
+  PetscScalar PETSC_UNUSED u = u0[0];
+  PetscScalar PETSC_UNUSED v = u0[1];
+  PetscScalar PETSC_UNUSED w = u0[2];
+  PetscScalar PETSC_UNUSED r = u0[3];
+
+  PetscScalar PETSC_UNUSED u_x = u1[0][0], u_y = u1[0][1];
+  PetscScalar PETSC_UNUSED v_x = u1[1][0], v_y = u1[1][1];
+  PetscScalar PETSC_UNUSED w_x = u1[2][0], w_y = u1[2][1];
+  PetscScalar PETSC_UNUSED r_x = u1[3][0], r_y = u1[3][1];
+
+  PetscScalar PETSC_UNUSED u_xx = u2[0][0][0], u_yy = u2[0][1][1];
+  PetscScalar PETSC_UNUSED v_xx = u2[1][0][0], v_yy = u2[1][1][1];
+  PetscScalar PETSC_UNUSED w_xx = u2[2][0][0], w_yy = u2[2][1][1];
+  PetscScalar PETSC_UNUSED r_xx = u2[3][0][0], r_yy = u2[3][1][1];
+
+  Fe[0] = u - Peaks(x,y);
+  Fe[1] = -(v_xx + v_yy) - 1.0;
+  Fe[2] = -(w_xx + w_yy) + w - 1.0;
+  Fe[3] = -(r_xx + r_yy) - 1*exp(r);
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "JacobianCollocation"
+PetscErrorCode JacobianCollocation(IGAPoint p,const PetscScalar *Ue,PetscScalar *Je,void *ctx)
+{
+  PetscInt nen = p->nen;
+
+  PetscScalar u0[4];
+  IGAPointFormValue(p,Ue,&u0[0]);
+  PetscScalar r = u0[3];
+
+  PetscReal N0[nen],N2[nen][2][2];
+  IGAPointFormShapeFuns(p,0,&N0[0]);
+  IGAPointFormShapeFuns(p,2,&N2[0][0][0]);
+
+  PetscInt a;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_xx = N2[a][0][0];
+    PetscReal Na_yy = N2[a][1][1];
+    Je[0*nen*4+a*4+0] = Na;
+    Je[1*nen*4+a*4+1] = -(Na_xx + Na_yy);
+    Je[2*nen*4+a*4+2] = -(Na_xx + Na_yy) + Na;
+    Je[3*nen*4+a*4+3] = -(Na_xx + Na_yy) - Na * 1*exp(r);
+  }
+  return 0;
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt N=16, p=2, C=PETSC_DECIDE;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","MultiComp2D Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-N","number of elements (along one dimension)",__FILE__,N,&N,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-p","polynomial order",__FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  if (C == PETSC_DECIDE) C = p-1;
+  ierr = IGAOptionsAlias("-N","16","-iga_elements");CHKERRQ(ierr);
+  ierr = IGAOptionsAlias("-p", "2","-iga_degree");CHKERRQ(ierr);
+  ierr = IGAOptionsAlias("-C","-1","-iga_continuity");CHKERRQ(ierr);
+  ierr = IGAOptionsAlias("-L","-1,+1","-iga_limits");CHKERRQ(ierr);
 
 #if PETSC_VERSION_(3,2,0)
   {
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,4);CHKERRQ(ierr);
 
-  IGAAxis axis0;
-  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisSetDegree(axis0,p);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,N,-1.0,1.0,C);CHKERRQ(ierr);
-  IGAAxis axis1;
-  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
-  ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
+  ierr = IGASetFieldName(iga,0,"L2Projection");CHKERRQ(ierr);
+  ierr = IGASetFieldName(iga,1,"Poisson");CHKERRQ(ierr);
+  ierr = IGASetFieldName(iga,2,"ReactionDiffusion");CHKERRQ(ierr);
+  ierr = IGASetFieldName(iga,3,"Bratu");CHKERRQ(ierr);
 
   IGABoundary bnd;
   PetscInt dir,side;
   for (dir=0; dir<2; dir++) {
     for (side=0; side<2; side++) {
+      PetscScalar field = 1;
       PetscScalar value = 1.0;
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,1,value);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
     }
   }
-  for (side=0, dir=0; dir<2; dir++) {
+  for (dir=0; dir<2; dir++) {
+    PetscScalar field = 2;
     PetscScalar value = 0.0;
-    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,2,value);CHKERRQ(ierr);
+    /*PetscScalar load  = 0.0;*/
+    ierr = IGAGetBoundary(iga,dir,side=0,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,dir,side=1,&bnd);CHKERRQ(ierr);
+    /*ierr = IGABoundarySetLoad(bnd,field,load);CHKERRQ(ierr);*/
+    ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
   }
   for (dir=0; dir<2; dir++) {
     for (side=0; side<2; side++) {
+      PetscScalar field = 3;
       PetscScalar value = 0.0;
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,3,value);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
     }
   }
 
-  ierr = IGASetUserFunction(iga,Function,0);CHKERRQ(ierr);
-  ierr = IGASetUserJacobian(iga,Jacobian,0);CHKERRQ(ierr);
-
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
+  if (!iga->collocation) {
+    ierr = IGASetUserFunction(iga,Function,0);CHKERRQ(ierr);
+    ierr = IGASetUserJacobian(iga,Jacobian,0);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserFunction(iga,FunctionCollocation,0);CHKERRQ(ierr);
+    ierr = IGASetUserJacobian(iga,JacobianCollocation,0);CHKERRQ(ierr);
+  }
+
   Vec x;
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   SNES snes;
   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
 
-  PetscBool flag = PETSC_FALSE;
-  PetscReal secs = -1;
-  ierr = PetscOptionsHasName(0,"-sleep",&flag);CHKERRQ(ierr);
-  ierr = PetscOptionsGetReal(0,"-sleep",&secs,0);CHKERRQ(ierr);
-  if (flag) {ierr = PetscSleep(secs);CHKERRQ(ierr);}
-
   ierr = PetscFinalize();CHKERRQ(ierr);
   return 0;
 }

File test/makefile

 	-@${MPIEXEC} -n 4 ./Test_SNES_2D ${OPTS} -N 16 -p 1 -ksp_type cg -iga_mat_type is -is_mat_type baij
 	-@${MPIEXEC} -n 4 ./Test_SNES_2D ${OPTS} -N 16 -p 1 -ksp_type cg -iga_mat_type is -is_mat_type sbaij
 	-@${MPIEXEC} -n 4 ./Test_SNES_2D ${OPTS} -N  8 -p 1 -ksp_type cg -iga_mat_type is -is_mat_type dense
+runex0e_1:
+	-@${MPIEXEC} -n 1 ./Test_SNES_2D ${OPTS} -p 2 -iga_collocation
+runex0e_4:
+	-@${MPIEXEC} -n 4 ./Test_SNES_2D ${OPTS} -p 4 -iga_collocation
 
 Test_SNES_2D = Test_SNES_2D.PETSc  \
 	       runex0a_1 runex0a_4 \
 	       runex0b_1 runex0b_4 \
 	       runex0c_1 runex0c_4 \
 	       runex0d_1 runex0d_4 \
+	       runex0e_1 runex0e_4 \
 	       Test_SNES_2D.rm
 
 
 	      ${MPIEXEC} -n 2 ./Test_SNES_2D -nox -malloc_debug -malloc_dump; \
 	      ${MPIEXEC} -n 3 ./Test_SNES_2D -nox -malloc_debug -malloc_dump; \
 	      ${MPIEXEC} -n 4 ./Test_SNES_2D -nox -malloc_debug -malloc_dump; \
+	      ${MPIEXEC} -n 9 ./Test_SNES_2D -nox -malloc_debug -malloc_dump; \
 	    fi; \
 	    ${OMAKE} Test_SNES_2D.rm PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}; \
 	   fi