Commits

Lisandro Dalcin committed acbe5e1

Rework IGAElement/IGAPoint preparing for changes on collocation

Comments (0)

Files changed (11)

demo/ClassicalShell.c

   PetscReal t  = user->t;
 
   // get geometry
-  PetscScalar *geom = p->parent->geometryX;
+  PetscScalar *geom = p->geometry;
   PetscScalar grad_g[3][2],hess_g[3][2][2];
   IGAShellInterpolate(p,1,3,geom,&grad_g[0][0]);
   IGAShellInterpolate(p,2,3,geom,&hess_g[0][0][0]);
   PetscInt dim;
   PetscInt nsd;
 
-  PetscInt  *mapping;   /*   [nen]      */
+  PetscInt  *mapping;  /*   [nen]      */
 
   PetscBool geometry;
   PetscBool rational;
-  PetscReal *geometryX; /*   [nen][nsd] */
-  PetscReal *geometryW; /*   [nen]      */
+  PetscReal *geometryX;/*   [nen][nsd] */
+  PetscReal *geometryW;/*   [nen]      */
 
   PetscReal *weight;   /*   [nqp]                     */
   PetscReal *detJac;   /*   [nqp]                     */
   IGA      parent;
   IGAPoint iterator;
 
+  PetscInt     nrows,*irows;
+  PetscInt     ncols,*icols;
+
   PetscInt     nfix;
   PetscInt    *ifix;
   PetscScalar *vfix;
   PetscScalar *xfix;
 
+  PetscInt    nval;
+  PetscScalar *wval[8];
   PetscInt    nvec;
   PetscScalar *wvec[8];
   PetscInt    nmat;
 PETSC_EXTERN PetscErrorCode IGAElementReset(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementInit(IGAElement element,IGA iga);
 
-PETSC_EXTERN PetscErrorCode IGAElementBegin(IGAElement element);
-PETSC_EXTERN PetscBool      IGAElementNext(IGAElement element);
-PETSC_EXTERN PetscErrorCode IGAElementEnd(IGAElement element);
+PETSC_EXTERN PetscErrorCode IGABeginElement(IGA iga,IGAElement *element);
+PETSC_EXTERN PetscBool      IGANextElement(IGA iga,IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAEndElement(IGA iga,IGAElement *element);
 
-PETSC_EXTERN PetscErrorCode IGAElementBuildFix(IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAElementInitPoint(IGAElement element,IGAPoint point);
+PETSC_EXTERN PetscErrorCode IGAElementBeginPoint(IGAElement element,IGAPoint *point);
+PETSC_EXTERN PetscBool      IGAElementNextPoint(IGAElement element,IGAPoint point);
+PETSC_EXTERN PetscErrorCode IGAElementEndPoint(IGAElement element,IGAPoint *point);
+
 PETSC_EXTERN PetscErrorCode IGAElementBuildMapping(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildGeometry(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildQuadrature(IGAElement element);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
 
+PETSC_EXTERN PetscErrorCode IGAElementGetWorkVal(IGAElement element,PetscScalar *U[]);
 PETSC_EXTERN PetscErrorCode IGAElementGetWorkVec(IGAElement element,PetscScalar *V[]);
 PETSC_EXTERN PetscErrorCode IGAElementGetWorkMat(IGAElement element,PetscScalar *M[]);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetValues(IGAElement element,const PetscScalar U[],PetscScalar u[]);
 
+PETSC_EXTERN PetscErrorCode IGAElementBuildFix(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementFixValues(IGAElement element,PetscScalar U[]);
 PETSC_EXTERN PetscErrorCode IGAElementFixFunction(IGAElement element,PetscScalar F[]);
 PETSC_EXTERN PetscErrorCode IGAElementFixJacobian(IGAElement element,PetscScalar J[]);
 
   PetscReal *point;    /*   [dim] */
   PetscReal *scale;    /*   [dim] */
-  
   PetscReal *basis[4]; /*0: [nen] */
                        /*1: [nen][dim] */
                        /*2: [nen][dim][dim] */
                        /*3: [nen][dim][dim][dim] */
 
-  PetscReal *detX;     /*   [nqp] */
+  PetscReal *geometry; /*   [nen][nsd] */
+  PetscReal *detX;     /*   [1] */
   PetscReal *gradX[2]; /*0: [nsd][dim] */
                        /*1: [dim][nsd] */
   PetscReal *shape[4]; /*0: [nen]  */
                        /*2: [nen][nsd][nsd] */
                        /*3: [nen][nsd][nsd][nsd] */
 
-  IGAElement  parent;
-
   PetscInt    nvec;
   PetscScalar *wvec[8];
   PetscInt    nmat;
 PETSC_EXTERN PetscErrorCode IGAPointCreate(IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAPointDestroy(IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAPointReset(IGAPoint point);
-PETSC_EXTERN PetscErrorCode IGAPointInit(IGAPoint point,IGAElement element);
-
-PETSC_EXTERN PetscErrorCode IGAPointBegin(IGAPoint point);
-PETSC_EXTERN PetscBool      IGAPointNext(IGAPoint point);
 
 PETSC_EXTERN PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index);
 PETSC_EXTERN PetscErrorCode IGAPointGetCount(IGAPoint point,PetscInt *count);
       ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);
     }
     ierr = IGABasisInit(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
+  }
+  ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
+
+  for (i=0; i<3; i++) {
     ierr = IGAColBasisInit(iga->colbasis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
   }
-  ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
   ierr = IGAColPointInit(iga->point_iterator,iga);CHKERRQ(ierr);
 
 
 PetscErrorCode IGAElementFreeWork(IGAElement element)
 {
   size_t i;
+  size_t MAX_WORK_VAL;
   size_t MAX_WORK_VEC;
   size_t MAX_WORK_MAT;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
+  MAX_WORK_VAL = sizeof(element->wval)/sizeof(PetscScalar*);
+  for (i=0; i<MAX_WORK_VAL; i++)
+    {ierr = PetscFree(element->wval[i]);CHKERRQ(ierr);}
+  element->nval = 0;
   MAX_WORK_VEC = sizeof(element->wvec)/sizeof(PetscScalar*);
   for (i=0; i<MAX_WORK_VEC; i++)
     {ierr = PetscFree(element->wvec[i]);CHKERRQ(ierr);}
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAElementFreeFix"
+static
+PetscErrorCode IGAElementFreeFix(IGAElement element)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  element->nfix = 0;
+  ierr = PetscFree(element->ifix);CHKERRQ(ierr);
+  ierr = PetscFree(element->vfix);CHKERRQ(ierr);
+  ierr = PetscFree(element->xfix);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAElementReset"
 PetscErrorCode IGAElementReset(IGAElement element)
 {
   PetscFunctionBegin;
   if (!element) PetscFunctionReturn(0);
   PetscValidPointer(element,1);
+  element->count =  0;
   element->index = -1;
 
   ierr = PetscFree(element->mapping);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[2]);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[3]);CHKERRQ(ierr);
 
+  element->nrows = 0; element->irows = PETSC_NULL;
+  element->ncols = 0; element->icols = PETSC_NULL;
+  ierr = IGAElementFreeFix(element);CHKERRQ(ierr);
   ierr = IGAElementFreeWork(element);CHKERRQ(ierr);
-  ierr = PetscFree(element->ifix);CHKERRQ(ierr);
-  ierr = PetscFree(element->vfix);CHKERRQ(ierr);
-  ierr = PetscFree(element->xfix);CHKERRQ(ierr);
   ierr = IGAPointReset(element->iterator);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
-  if (PetscUnlikely(!iga->setup)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetUp() first");
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,2);
+  IGACheckSetUp(iga,1);
   ierr = IGAElementReset(element);CHKERRQ(ierr);
   element->parent = iga;
 
     ierr = PetscMemzero(element->shape[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
   }
-  ierr = IGAPointInit(element->iterator,element);CHKERRQ(ierr);
+  { /* */
+    element->nrows = element->nen; element->irows = element->mapping;
+    element->ncols = element->nen; element->icols = element->mapping;
+  }
   { /* */
     PetscInt nen = element->nen;
     PetscInt dof = element->dof;
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->vfix);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->xfix);CHKERRQ(ierr);
   }
+  ierr = IGAElementInitPoint(element,element->iterator);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementBegin"
-PetscErrorCode IGAElementBegin(IGAElement element)
+#define __FUNCT__ "IGABeginElement"
+PetscErrorCode IGABeginElement(IGA iga,IGAElement *element)
 {
-  IGA            iga;
   PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  iga = element->parent;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
-  if (PetscUnlikely(!iga->setup)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetUp() first");
-  element->index = -1;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(element,2);
+  IGACheckSetUp(iga,1);
+  *element = iga->iterator;
+  (*element)->index = -1;
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementNext"
-PetscBool IGAElementNext(IGAElement element)
+#define __FUNCT__ "IGANextElement"
+PetscBool IGANextElement(IGA iga,IGAElement element)
 {
   PetscInt i,dim  = element->dim;
   PetscInt *start = element->start;
   PetscInt index,coord;
   PetscErrorCode ierr;
   /* */
+  element->nval = 0;
   element->nvec = 0;
   element->nmat = 0;
   /* */
     index = (index - coord) / width[i];
     ID[i] = coord + start[i];
   }
-#undef  CHKERR_RETURN
-#define CHKERR_RETURN(n,r) do{if(PetscUnlikely(n)){CHKERRCONTINUE(n);return(r);}}while(0)
-  ierr = IGAElementBuildMapping(element);  CHKERR_RETURN(ierr,PETSC_FALSE);
-  ierr = IGAElementBuildGeometry(element); CHKERR_RETURN(ierr,PETSC_FALSE);
-  ierr = IGAElementBuildFix(element);      CHKERR_RETURN(ierr,PETSC_FALSE);
-#undef  CHKERR_RETURN_FALSE
+#undef  CHKERRRETURN
+#define CHKERRRETURN(n,r) do{if(PetscUnlikely(n)){CHKERRCONTINUE(n);return(r);}}while(0)
+  ierr = IGAElementBuildMapping(element);  CHKERRRETURN(ierr,PETSC_FALSE);
+  ierr = IGAElementBuildGeometry(element); CHKERRRETURN(ierr,PETSC_FALSE);
+  ierr = IGAElementBuildFix(element);      CHKERRRETURN(ierr,PETSC_FALSE);
+#undef  CHKERRRETURN
   return PETSC_TRUE;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementEnd"
-PetscErrorCode IGAElementEnd(IGAElement element)
+#define __FUNCT__ "IGAEndElement"
+PetscErrorCode IGAEndElement(IGA iga,IGAElement *element)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(element,2);
+  PetscValidPointer(*element,2);
+  if (PetscUnlikely((*element)->index != -1)) {
+    (*element)->index = -1;
+    *element = PETSC_NULL;
+    PetscFunctionReturn(PETSC_ERR_PLIB);
+  }
+  *element = PETSC_NULL;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementInitPoint"
+PetscErrorCode IGAElementInitPoint(IGAElement element,IGAPoint point)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  PetscValidPointer(point,2);
+  ierr = IGAPointReset(point);CHKERRQ(ierr);
+
+  point->nen = element->nen;
+  point->dof = element->dof;
+  point->dim = element->dim;
+  point->nsd = element->nsd;
+
+  point->count = element->nqp;
+  point->index = -1;
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementBeginPoint"
+PetscErrorCode IGAElementBeginPoint(IGAElement element,IGAPoint *point)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  PetscValidPointer(point,2);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
+  ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
+  *point = element->iterator;
+  (*point)->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementNextPoint"
+PetscBool IGAElementNextPoint(IGAElement element,IGAPoint point)
+{
+  PetscInt nen = point->nen;
+  PetscInt dim = point->dim;
+  PetscInt nsd = point->nsd;
+  PetscInt index;
+  /*PetscValidPointer(element,1);*/
+  /*PetscValidPointer(point,2);*/
+  /* */
+  point->nvec = 0;
+  point->nmat = 0;
+  /* */
+  index = ++point->index;
+  if (PetscUnlikely(index == 0))            goto start;
+  if (PetscUnlikely(index >= point->count)) goto stop;
+
+  point->weight   += 1;
+  point->detJac   += 1;
+
+  point->point    += dim;
+  point->scale    += dim;
+  point->basis[0] += nen;
+  point->basis[1] += nen*dim;
+  point->basis[2] += nen*dim*dim;
+  point->basis[3] += nen*dim*dim*dim;
+
+  point->detX     += 1;
+  point->gradX[0] += dim*dim;
+  point->gradX[1] += dim*dim;
+  point->shape[0] += nen;
+  point->shape[1] += nen*dim;
+  point->shape[2] += nen*dim*dim;
+  point->shape[3] += nen*dim*dim*dim;
+
+  return PETSC_TRUE;
+
+ start:
+
+  point->weight   = element->weight;
+  point->detJac   = element->detJac;
+
+  point->point    = element->point;
+  point->scale    = element->scale;
+  point->basis[0] = element->basis[0];
+  point->basis[1] = element->basis[1];
+  point->basis[2] = element->basis[2];
+  point->basis[3] = element->basis[3];
+
+  if (element->geometry)
+    point->geometry = element->geometryX;
+  else 
+    point->geometry = PETSC_NULL;
+  if (element->geometry && dim == nsd) { /* XXX */
+    point->detX     = element->detX;
+    point->gradX[0] = element->gradX[0];
+    point->gradX[1] = element->gradX[1];
+    point->shape[0] = element->shape[0];
+    point->shape[1] = element->shape[1];
+    point->shape[2] = element->shape[2];
+    point->shape[3] = element->shape[3];
+  } else {
+    point->shape[0] = element->basis[0];
+    point->shape[1] = element->basis[1];
+    point->shape[2] = element->basis[2];
+    point->shape[3] = element->basis[3];
+  }
+
+  return PETSC_TRUE;
+
+ stop:
+
+  point->index = -1;
+  return PETSC_FALSE;
+
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementEndPoint"
+PetscErrorCode IGAElementEndPoint(IGAElement element,IGAPoint *point)
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
-  if (PetscUnlikely(element->index != -1)) {
-    element->index = -1;
+  PetscValidPointer(point,2);
+  PetscValidPointer(*point,2);
+  if (PetscUnlikely((*point)->index != -1)) {
+    (*point)->index = -1;
     PetscFunctionReturn(PETSC_ERR_PLIB);
   }
+  *point = PETSC_NULL;
   PetscFunctionReturn(0);
 }
 
 
 /* -------------------------------------------------------------------------- */
 
-static void AddFix(IGAElement element,PetscInt dir,PetscInt side,PetscInt a)
-{
-  IGABoundary b = element->parent->boundary[dir][side];
-  PetscInt dof = element->dof;
-  PetscInt nfix = element->nfix;
-  PetscInt *ifix = element->ifix;
-  PetscScalar *vfix = element->vfix;
-  PetscInt j,k,n = b->nbc;
-  for (k=0; k<n; k++) {
-    PetscInt idx = a*dof + b->field[k];
-    PetscScalar val = b->value[k];
-    for (j=0; j<nfix; j++)
-      if (ifix[j] == idx) break;
-    if (j==nfix) nfix++;
-    ifix[j] = idx;
-    vfix[j] = val;
-  }
-  element->nfix = nfix;
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementBuildFix"
-PetscErrorCode IGAElementBuildFix(IGAElement element)
-{
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  if (PetscUnlikely(element->index < 0))
-    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  {
-    IGA      iga = element->parent;
-    IGABasis *BD = iga->basis;
-    PetscInt *ID = element->ID;
-    PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
-    PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
-    PetscBool onboundary = PETSC_FALSE;
-    PetscInt i,dim = element->dim;
-    for (i=0; i<dim; i++) {
-      PetscInt e = BD[i]->nel-1; /* last element */
-      PetscInt n = BD[i]->nnp-1; /* last node */
-      if (iga->axis[i]->periodic) continue; /* XXX */
-      if (ID[i] == 0) { A0[i] = 0; onboundary = PETSC_TRUE; }
-      if (ID[i] == e) { A1[i] = n; onboundary = PETSC_TRUE; }
-    }
-    element->nfix = 0;
-    if (onboundary) {
-      PetscInt ia, inen = BD[0]->nen, ioffset = BD[0]->offset[ID[0]];
-      PetscInt ja, jnen = BD[1]->nen, joffset = BD[1]->offset[ID[1]];
-      PetscInt ka, knen = BD[2]->nen, koffset = BD[2]->offset[ID[2]];
-      PetscInt a = 0;
-      for (ka=0; ka<knen; ka++) {
-        for (ja=0; ja<jnen; ja++) {
-          for (ia=0; ia<inen; ia++) {
-            PetscInt iA = ioffset + ia;
-            PetscInt jA = joffset + ja;
-            PetscInt kA = koffset + ka;
-            /**/ if (iA == A0[0]) AddFix(element,0,0,a);
-            else if (iA == A1[0]) AddFix(element,0,1,a);
-            /**/ if (jA == A0[1]) AddFix(element,1,0,a);
-            else if (jA == A1[1]) AddFix(element,1,1,a);
-            /**/ if (kA == A0[2]) AddFix(element,2,0,a);
-            else if (kA == A1[2]) AddFix(element,2,1,a);
-            a++;
-          }
-        }
-      }
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementBuildMapping"
 PetscErrorCode IGAElementBuildMapping(IGAElement element)
 /* -------------------------------------------------------------------------- */
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAElementGetWorkVal"
+PetscErrorCode IGAElementGetWorkVal(IGAElement element,PetscScalar *U[])
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  PetscValidPointer(U,2);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  {
+    size_t MAX_WORK_VAL = sizeof(element->wval)/sizeof(PetscScalar*);
+    PetscInt m = element->nen * element->dof;
+    if (PetscUnlikely(element->nval >= (PetscInt)MAX_WORK_VAL))
+      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work values requested");
+    if (PetscUnlikely(!element->wval[element->nval])) {
+      ierr = PetscMalloc1(m,PetscScalar,&element->wval[element->nval]);CHKERRQ(ierr);
+    }
+    *U = element->wval[element->nval++];
+    ierr = PetscMemzero(*U,m*sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetWorkVec"
 PetscErrorCode IGAElementGetWorkVec(IGAElement element,PetscScalar *V[])
 {
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
     size_t MAX_WORK_VEC = sizeof(element->wvec)/sizeof(PetscScalar*);
-    PetscInt n = element->nen * element->dof;
+    PetscInt m = element->nrows * element->dof;
     if (PetscUnlikely(element->nvec >= (PetscInt)MAX_WORK_VEC))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work vectors requested");
     if (PetscUnlikely(!element->wvec[element->nvec])) {
-      ierr = PetscMalloc1(n,PetscScalar,&element->wvec[element->nvec]);CHKERRQ(ierr);
+      ierr = PetscMalloc1(m,PetscScalar,&element->wvec[element->nvec]);CHKERRQ(ierr);
     }
     *V = element->wvec[element->nvec++];
-    ierr = PetscMemzero(*V,n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = PetscMemzero(*V,m*sizeof(PetscScalar));CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
     size_t MAX_WORK_MAT = sizeof(element->wmat)/sizeof(PetscScalar*);
-    PetscInt n = element->nen * element->dof;
+    PetscInt m = element->nrows * element->dof;
+    PetscInt n = element->ncols * element->dof;
     if (PetscUnlikely(element->nmat >= (PetscInt)MAX_WORK_MAT))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work matrices requested");
     if (PetscUnlikely(!element->wmat[element->nmat])) {
-      ierr = PetscMalloc1(n*n,PetscScalar,&element->wmat[element->nmat]);CHKERRQ(ierr);
+      ierr = PetscMalloc1(m*n,PetscScalar,&element->wmat[element->nmat]);CHKERRQ(ierr);
     }
     *M = element->wmat[element->nmat++];
-    ierr = PetscMemzero(*M,n*n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = PetscMemzero(*M,m*n*sizeof(PetscScalar));CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetValues"
-PetscErrorCode IGAElementGetValues(IGAElement element,const PetscScalar arrayU[], PetscScalar U[])
+PetscErrorCode IGAElementGetValues(IGAElement element,const PetscScalar arrayU[],PetscScalar U[])
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscFunctionReturn(0);
 }
 
+static void AddFix(IGAElement element,PetscInt dir,PetscInt side,PetscInt a)
+{
+  IGABoundary b = element->parent->boundary[dir][side];
+  PetscInt dof = element->dof;
+  PetscInt nfix = element->nfix;
+  PetscInt *ifix = element->ifix;
+  PetscScalar *vfix = element->vfix;
+  PetscInt j,k,n = b->nbc;
+  for (k=0; k<n; k++) {
+    PetscInt idx = a*dof + b->field[k];
+    PetscScalar val = b->value[k];
+    for (j=0; j<nfix; j++)
+      if (ifix[j] == idx) break;
+    if (j==nfix) nfix++;
+    ifix[j] = idx;
+    vfix[j] = val;
+  }
+  element->nfix = nfix;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementBuildFix"
+PetscErrorCode IGAElementBuildFix(IGAElement element)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  {
+    IGA      iga = element->parent;
+    IGABasis *BD = iga->basis;
+    PetscInt *ID = element->ID;
+    PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
+    PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
+    PetscBool onboundary = PETSC_FALSE;
+    PetscInt i,dim = element->dim;
+    for (i=0; i<dim; i++) {
+      PetscInt e = BD[i]->nel-1; /* last element */
+      PetscInt n = BD[i]->nnp-1; /* last node    */
+      if (iga->axis[i]->periodic) continue; /* XXX */
+      if (ID[i] == 0) { A0[i] = 0; onboundary = PETSC_TRUE; }
+      if (ID[i] == e) { A1[i] = n; onboundary = PETSC_TRUE; }
+    }
+    element->nfix = 0;
+    if (onboundary) {
+      PetscInt ia, inen = BD[0]->nen, ioffset = BD[0]->offset[ID[0]];
+      PetscInt ja, jnen = BD[1]->nen, joffset = BD[1]->offset[ID[1]];
+      PetscInt ka, knen = BD[2]->nen, koffset = BD[2]->offset[ID[2]];
+      PetscInt a = 0;
+      for (ka=0; ka<knen; ka++) {
+        for (ja=0; ja<jnen; ja++) {
+          for (ia=0; ia<inen; ia++) {
+            PetscInt iA = ioffset + ia;
+            PetscInt jA = joffset + ja;
+            PetscInt kA = koffset + ka;
+            /**/ if (iA == A0[0]) AddFix(element,0,0,a);
+            else if (iA == A1[0]) AddFix(element,0,1,a);
+            /**/ if (jA == A0[1]) AddFix(element,1,0,a);
+            else if (jA == A1[1]) AddFix(element,1,1,a);
+            /**/ if (kA == A0[2]) AddFix(element,2,0,a);
+            else if (kA == A1[2]) AddFix(element,2,1,a);
+            a++;
+          }
+        }
+      }
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementFixValues"
 PetscErrorCode IGAElementFixValues(IGAElement element,PetscScalar U[])
   PetscValidPointer(element,1);
   PetscValidScalarPointer(U,2);
   {
-    PetscInt k,nfix=element->nfix;
-    for (k=0; k<nfix; k++) {
-      PetscInt i = element->ifix[k];
-      element->xfix[k] = U[i];
-      U[i] = element->vfix[k];
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      PetscInt k = element->ifix[f];
+      element->xfix[f] = U[k];
+      U[k] = element->vfix[f];
     }
   }
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAElementFixSystem"
+PetscErrorCode IGAElementFixSystem(IGAElement element,PetscScalar K[],PetscScalar F[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  PetscValidScalarPointer(K,2);
+  PetscValidScalarPointer(F,3);
+  {
+    PetscInt M = element->nrows*element->dof;
+    PetscInt N = element->ncols*element->dof;
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      PetscInt    k = element->ifix[f];
+      PetscScalar v = element->vfix[f];
+      PetscInt i,j;
+      for (i=0; i<M; i++)
+        F[i] -= K[i*N+k] * v;
+      for (i=0; i<M; i++) 
+        K[i*N+k] = 0.0;
+      for (j=0; j<N; j++) 
+        K[k*N+j] = 0.0;
+      K[k*N+k] = 1.0;
+      F[k] = v;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+#undef  __FUNCT__
 #define __FUNCT__ "IGAElementFixFunction"
 PetscErrorCode IGAElementFixFunction(IGAElement element,PetscScalar F[])
 {
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
   {
-    PetscInt k,nfix=element->nfix;
-    for (k=0; k<nfix; k++) {
-      PetscInt i = element->ifix[k];
-      PetscScalar u = element->xfix[k];
-      PetscScalar u0 = element->vfix[k];
-      F[i] = u - u0;
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      PetscInt    k  = element->ifix[f];
+      PetscScalar u  = element->xfix[f];
+      PetscScalar u0 = element->vfix[f];
+      F[k] = u - u0;
     }
   }
   PetscFunctionReturn(0);
   PetscValidPointer(element,1);
   PetscValidScalarPointer(J,2);
   {
-    PetscInt N = element->nen*element->dof;
-    PetscInt k,nfix=element->nfix;
-    for (k=0; k<nfix; k++) {
-      PetscInt j,i = element->ifix[k];
-      for (j=0; j<N; j++)
-        J[i*N+j] = J[j*N+i] = 0.0;
-      J[i*N+i] = 1.0;
-    }
-  }
-  PetscFunctionReturn(0);
-}
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementFixSystem"
-PetscErrorCode IGAElementFixSystem(IGAElement element,PetscScalar K[],PetscScalar F[])
-{
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  PetscValidScalarPointer(K,2);
-  PetscValidScalarPointer(F,3);
-  {
-    PetscInt N = element->nen*element->dof;
-    PetscInt k,nfix=element->nfix;
-    for (k=0; k<nfix; k++) {
-      PetscInt j,i = element->ifix[k];
-      PetscScalar v = element->vfix[k];
-      for (j=0; j<N; j++) {
-        F[j] -= K[j*N+i] * v;
-      }
-      for (j=0; j<N; j++)
-        K[i*N+j] = K[j*N+i] = 0.0;
-      K[i*N+i] = 1.0;
-      F[i] = v;
+    PetscInt M = element->nrows*element->dof;
+    PetscInt N = element->ncols*element->dof;
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      PetscInt k = element->ifix[f];
+      PetscInt i,j;
+      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);
 #define __FUNCT__ "IGAElementAssembleVec"
 PetscErrorCode IGAElementAssembleVec(IGAElement element,const PetscScalar F[],Vec vec)
 {
-  PetscInt       nn,*ii;
+  PetscInt       mm,*ii;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
   PetscValidHeaderSpecific(vec,VEC_CLASSID,3);
-  nn = element->nen;
-  ii = element->mapping;
+  mm = element->nrows; ii = element->irows;
   if (element->dof == 1) {
-    ierr = VecSetValuesLocal(vec,nn,ii,F,ADD_VALUES);CHKERRQ(ierr);
+    ierr = VecSetValuesLocal(vec,mm,ii,F,ADD_VALUES);CHKERRQ(ierr);
   } else {
-    ierr = VecSetValuesBlockedLocal(vec,nn,ii,F,ADD_VALUES);CHKERRQ(ierr);
+    ierr = VecSetValuesBlockedLocal(vec,mm,ii,F,ADD_VALUES);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 #define __FUNCT__ "IGAElementAssembleMat"
 PetscErrorCode IGAElementAssembleMat(IGAElement element,const PetscScalar K[],Mat mat)
 {
-  PetscInt       nn,*ii;
+  PetscInt       mm,*ii;
+  PetscInt       nn,*jj;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(K,2);
   PetscValidHeaderSpecific(mat,MAT_CLASSID,3);
-  nn = element->nen;
-  ii = element->mapping;
+  mm = element->nrows; ii = element->irows;
+  nn = element->ncols; jj = element->icols;
   if (element->dof == 1) {
-    ierr = MatSetValuesLocal(mat,nn,ii,nn,ii,K,ADD_VALUES);CHKERRQ(ierr);
+    ierr = MatSetValuesLocal(mat,mm,ii,nn,jj,K,ADD_VALUES);CHKERRQ(ierr);
   } else {
-    ierr = MatSetValuesBlockedLocal(mat,nn,ii,nn,ii,K,ADD_VALUES);CHKERRQ(ierr);
+    ierr = MatSetValuesBlockedLocal(mat,mm,ii,nn,jj,K,ADD_VALUES);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }

src/petigaftn.F90

      type(C_PTR) :: point
      type(C_PTR) :: scale
      type(C_PTR) :: basis(0:3)
+     type(C_PTR) :: geometry
      type(C_PTR) :: detX
      type(C_PTR) :: gradX(0:1)
      type(C_PTR) :: shape(0:3)
   ierr = MatZeroEntries(matA);CHKERRQ(ierr);
   ierr = VecZeroEntries(vecB);CHKERRQ(ierr);
 
+  /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *A, *B;
     ierr = IGAElementGetWorkMat(element,&A);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVec(element,&B);CHKERRQ(ierr);
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    /* Quadrature loop */
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *K, *F;
       ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
       ierr = IGAPointGetWorkVec(point,&F);CHKERRQ(ierr);
       ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
       ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
+    /* */
     ierr = IGAElementFixSystem(element,A,B);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,A,matA);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
 
   ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     lwork = (info==0) ? (PetscBLASInt)work[0] : m*128;
     ierr = PetscMalloc1(lwork,PetscScalar,&work);CHKERRQ(ierr);
 
-    ierr = IGAElementBegin(element);CHKERRQ(ierr);
-    while (IGAElementNext(element)) {
+    ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+    while (IGANextElement(iga,element)) {
       ierr = IGAElementGetMapping(element,&nen,&mapping);CHKERRQ(ierr);
       m = n = ComputeOwnedGlobalIndices(ltogmap,start,end,dof,nen,mapping,indices);
       /* get element matrix from global matrix */
       ierr = MatSetValues(B,n,indices,n,indices,values,ADD_VALUES);CHKERRQ(ierr);
       ierr = PetscLogFlops(n*n);CHKERRQ(ierr);
     }
-    ierr = IGAElementEnd(element);CHKERRQ(ierr);
+    ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
 
     ierr = ISLocalToGlobalMappingRestoreIndices(map,&ltogmap);CHKERRQ(ierr);
     ierr = PetscFree2(indices,values);CHKERRQ(ierr);

src/petigapoint.c

   PetscFunctionBegin;
   if (!point) PetscFunctionReturn(0);
   PetscValidPointer(point,1);
+  point->count =  0;
   point->index = -1;
   ierr = IGAPointFreeWork(point);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAPointInit"
-PetscErrorCode IGAPointInit(IGAPoint point,IGAElement element)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(point,1);
-  PetscValidPointer(element,0);
-  ierr = IGAPointReset(point);CHKERRQ(ierr);
-  point->parent = element;
-
-  point->nen = element->nen;
-  point->dof = element->dof;
-  point->dim = element->dim;
-  point->nsd = element->nsd;
-
-  point->count = element->nqp;
-  point->index = -1;
-
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAPointBegin"
-PetscErrorCode IGAPointBegin(IGAPoint point)
-{
-  IGAElement     element;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(point,1);
-  element = point->parent;
-  if (PetscUnlikely(element->index < 0))
-    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  point->index = -1;
-  ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
-  ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAPointNext"
-PetscBool IGAPointNext(IGAPoint point)
-{
-  IGAElement element;
-  PetscInt nen = point->nen;
-  PetscInt dim = point->dim;
-  PetscInt nsd = point->nsd;
-  PetscInt index;
-  /* */
-  point->nvec = 0;
-  point->nmat = 0;
-  /* */
-  index = ++point->index;
-  if (PetscUnlikely(index == 0))            goto start;
-  if (PetscUnlikely(index >= point->count)) goto stop;
-
-  point->weight   += 1;
-  point->detJac   += 1;
-
-  point->point    += dim;
-  point->scale    += dim;
-  point->basis[0] += nen;
-  point->basis[1] += nen*dim;
-  point->basis[2] += nen*dim*dim;
-  point->basis[3] += nen*dim*dim*dim;
-
-  point->detX     += 1;
-  point->gradX[0] += dim*dim;
-  point->gradX[1] += dim*dim;
-  point->shape[0] += nen;
-  point->shape[1] += nen*dim;
-  point->shape[2] += nen*dim*dim;
-  point->shape[3] += nen*dim*dim*dim;
-
-  return PETSC_TRUE;
-
- start:
-
-  element = point->parent;
-
-  point->weight   = element->weight;
-  point->detJac   = element->detJac;
-
-  point->point    = element->point;
-  point->scale    = element->scale;
-  point->basis[0] = element->basis[0];
-  point->basis[1] = element->basis[1];
-  point->basis[2] = element->basis[2];
-  point->basis[3] = element->basis[3];
-
-  if (element->geometry && dim == nsd) { /* XXX */
-    point->detX     = element->detX;
-    point->gradX[0] = element->gradX[0];
-    point->gradX[1] = element->gradX[1];
-    point->shape[0] = element->shape[0];
-    point->shape[1] = element->shape[1];
-    point->shape[2] = element->shape[2];
-    point->shape[3] = element->shape[3];
-  } else {
-    point->shape[0] = element->basis[0];
-    point->shape[1] = element->basis[1];
-    point->shape[2] = element->basis[2];
-    point->shape[3] = element->basis[3];
-  }
-  return PETSC_TRUE;
-
- stop:
-
-  point->index = -1;
-  return PETSC_FALSE;
-
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAPointGetParent"
-PetscErrorCode IGAPointGetParent(IGAPoint point,IGAElement *parent)
-{
-  PetscFunctionBegin;
-  PetscValidPointer(point,1);
-  PetscValidPointer(parent,2);
-  *parent = point->parent;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
 #define __FUNCT__ "IGAPointGetIndex"
 PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index)
 {
   PetscFunctionBegin;
   PetscValidPointer(p,1);
   PetscValidRealPointer(x,2);
-  if (p->parent->geometry) {
-    const PetscReal *X = p->parent->geometryX;
+  if (p->geometry) {
+    const PetscReal *X = p->geometry;
     IGA_GetGeomMap(p->nen,p->nsd,p->shape[0],X,x);
   } else {
     PetscInt i,dim = p->dim;
   PetscFunctionBegin;
   PetscValidPointer(p,1);
   PetscValidRealPointer(F,2);
-  if (p->parent->geometry) {
+  if (p->geometry) {
     PetscInt a,dim = p->dim;
     PetscInt i,nsd = p->nsd;
     const PetscReal *L = p->scale;
     if (dim == nsd) {
       (void)PetscMemcpy(F,p->gradX[0],nsd*dim*sizeof(PetscReal));
     } else {
-      const PetscReal *X = p->parent->geometryX;
+      const PetscReal *X = p->geometry;
       IGA_GetGradGeomMap(p->nen,nsd,dim,p->basis[1],X,F);
     }
     for (i=0; i<nsd; i++)
   PetscFunctionBegin;
   PetscValidPointer(p,1);
   PetscValidRealPointer(G,2);
-  if (p->parent->geometry) {
+  if (p->geometry) {
     PetscInt a,dim = p->dim;
     PetscInt i,nsd = p->nsd;
     const PetscReal *L = p->scale;
     if (dim == nsd) {
       (void)PetscMemcpy(G,p->gradX[1],dim*nsd*sizeof(PetscReal));
     } else {
-      const PetscReal *X = p->parent->geometryX;
+      const PetscReal *X = p->geometry;
       IGA_GetInvGradGeomMap(p->nen,nsd,dim,p->basis[1],X,G);
     }
     for (a=0; a<dim; a++)
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during point loop");
   {
     size_t MAX_WORK_VEC = sizeof(point->wvec)/sizeof(PetscScalar*);
-    PetscInt n = point->nen * point->dof;
+    PetscInt m = point->nen * point->dof;
     if (PetscUnlikely(point->nvec >= (PetscInt)MAX_WORK_VEC))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work vectors requested");
     if (PetscUnlikely(!point->wvec[point->nvec])) {
-      ierr = PetscMalloc1(n,PetscScalar,&point->wvec[point->nvec]);CHKERRQ(ierr);
+      ierr = PetscMalloc1(m,PetscScalar,&point->wvec[point->nvec]);CHKERRQ(ierr);
     }
     *V = point->wvec[point->nvec++];
-    ierr = PetscMemzero(*V,n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = PetscMemzero(*V,m*sizeof(PetscScalar));CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during point loop");
   {
     size_t MAX_WORK_MAT = sizeof(point->wmat)/sizeof(PetscScalar*);
+    PetscInt m = point->nen * point->dof;
     PetscInt n = point->nen * point->dof;
     if (PetscUnlikely(point->nmat >= (PetscInt)MAX_WORK_MAT))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work matrices requested");
     if (PetscUnlikely(!point->wmat[point->nmat])) {
-      ierr = PetscMalloc1(n*n,PetscScalar,&point->wmat[point->nmat]);CHKERRQ(ierr);
+      ierr = PetscMalloc1(m*n,PetscScalar,&point->wmat[point->nmat]);CHKERRQ(ierr);
     }
     *M = point->wmat[point->nmat++];
-    ierr = PetscMemzero(*M,n*n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = PetscMemzero(*M,m*n*sizeof(PetscScalar));CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 #define __FUNCT__ "IGAPointAddVec"
 PetscErrorCode IGAPointAddVec(IGAPoint point,const PetscScalar f[],PetscScalar F[])
 {
-  PetscInt       n;
+  PetscInt       m;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(point,1);
   PetscValidScalarPointer(f,2);
   PetscValidScalarPointer(F,3);
-  n = point->nen*point->dof;
-  ierr = IGAPointAddArray(point,n,f,F);CHKERRQ(ierr);
+  m = point->nen * point->dof;
+  ierr = IGAPointAddArray(point,m,f,F);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   PetscValidPointer(point,1);
   PetscValidScalarPointer(k,2);
   PetscValidScalarPointer(K,3);
-  m = n = point->nen*point->dof;
+  m = point->nen * point->dof;
+  n = point->nen * point->dof;
   ierr = IGAPointAddArray(point,m*n,k,K);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormScalar,iga,vecU,0,0);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *U;
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       ierr = PetscMemzero(workS,n*sizeof(PetscScalar));CHKERRQ(ierr);
       ierr = Scalar(point,U,n,workS,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddArray(point,n,workS,localS);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormScalar,iga,vecU,0,0);CHKERRQ(ierr);
 
   /* Restore local vector U and array */
 
   /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecU,vecF,0);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *U, *F;
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *R;
       ierr = IGAPointGetWorkVec(point,&R);CHKERRQ(ierr);
       ierr = Function(point,U,R,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddVec(point,R,F);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormFunction,iga,vecU,vecF,0);CHKERRQ(ierr);
 
   /* Restore local vector U and array */
 
   /* Element Loop */
   ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecU,matJ,0);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *U, *J;
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *K;
       ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
       ierr = Jacobian(point,U,K,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddMat(point,K,J);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormJacobian,iga,vecU,matJ,0);CHKERRQ(ierr);
 
   /* Restore local vector U and array */
 
   /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *V, *U, *F;
-    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *R;
       ierr = IGAPointGetWorkVec(point,&R);CHKERRQ(ierr);
       ierr = IFunction(point,dt,a,V,t,U,R,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddVec(point,R,F);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
 
   /* Restore local vectors V,U and arrays */
 
   /* Element Loop */
   ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *V, *U, *J;
-    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *K;
       ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
       ierr = IJacobian(point,dt,a,V,t,U,K,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddMat(point,K,J);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
 
   /* Get local vectors V,U and arrays */
 
   /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *V, *U, *U0, *F;
-    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U0);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U0);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU0,U0);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U0);CHKERRQ(ierr); /* XXX */
-    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);  /* XXX */
+    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *R;
       ierr = IGAPointGetWorkVec(point,&R);CHKERRQ(ierr);
       ierr = IEFunction(point,dt,a,V,t,U,t0,U0,R,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddVec(point,R,F);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
 
   /* Restore local vectors V,U,U0 and arrays */
 
   /* Element Loop */
   ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
-  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-  ierr = IGAElementBegin(element);CHKERRQ(ierr);
-  while (IGAElementNext(element)) {
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
     PetscScalar *V, *U, *U0, *J;
-    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&U0);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
-    /* */
+    ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(element,&U0);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU0,U0);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U0);CHKERRQ(ierr); /* XXX */
-    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);  /* XXX */
+    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
     /* Quadrature loop */
-    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
-    ierr = IGAPointBegin(point);CHKERRQ(ierr);
-    while (IGAPointNext(point)) {
+    ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+    while (IGAElementNextPoint(element,point)) {
       PetscScalar *K;
       ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
       ierr = IEJacobian(point,dt,a,V,t,U,t0,U0,K,ctx);CHKERRQ(ierr);
       ierr = IGAPointAddMat(point,K,J);CHKERRQ(ierr);
     }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
     /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
   }
-  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
   ierr = PetscLogEventEnd(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
 
   /* Restore local vectors V,U,U0 and arrays */