Commits

Lisandro Dalcin committed 8f8bb27

Rework and complete implementation of boundary integrals

Comments (0)

Files changed (8)

   PetscReal *point;   /* [nel][nqp]           */
   PetscReal *value;   /* [nel][nqp][nen][d+1] */
 
+  PetscInt   bnd_offset[2];
   PetscReal  bnd_detJ[2];
   PetscReal  bnd_weight[2];
   PetscReal  bnd_point[2];
 struct _IGAUserOps {
   /**/
   IGAUserSystem     System;
-  void              *SysCtx;  
+  void              *SysCtx;
   /**/
   IGAUserFunction   Function;
   void              *FunCtx;
 PETSC_EXTERN PetscErrorCode IGABoundaryClear(IGABoundary boundary);
 PETSC_EXTERN PetscErrorCode IGABoundarySetValue(IGABoundary boundary,PetscInt field,PetscScalar value);
 PETSC_EXTERN PetscErrorCode IGABoundarySetLoad (IGABoundary boundary,PetscInt field,PetscScalar value);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserSystem(IGABoundary boundary,IGAUserSystem System,void *SysCtx);
+
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserSystem    (IGABoundary boundary,IGAUserSystem     System,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserFunction  (IGABoundary boundary,IGAUserFunction   Function,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserJacobian  (IGABoundary boundary,IGAUserJacobian   Jacobian,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIFunction (IGABoundary boundary,IGAUserIFunction  IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIJacobian (IGABoundary boundary,IGAUserIJacobian  IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIFunction2(IGABoundary boundary,IGAUserIFunction2 IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIJacobian2(IGABoundary boundary,IGAUserIJacobian2 IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIEFunction(IGABoundary boundary,IGAUserIEFunction IEFunction,void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserIEJacobian(IGABoundary boundary,IGAUserIEJacobian IEJacobian,void *ctx);
+
 
 typedef struct _p_IGA *IGA;
 
 PETSC_EXTERN PetscErrorCode IGAGetLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 PETSC_EXTERN PetscErrorCode IGARestoreLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 
-PETSC_EXTERN PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *SysCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *FunCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  void *JacCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIFunction (IGA iga,IGAUserIFunction  IFunction, void *FunCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIJacobian (IGA iga,IGAUserIJacobian  IJacobian, void *JacCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIFunction2(IGA iga,IGAUserIFunction2 IFunction, void *FunCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIJacobian2(IGA iga,IGAUserIJacobian2 IJacobian, void *JacCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *FunCtx);
-PETSC_EXTERN PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *JacCtx);
+PETSC_EXTERN PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIFunction (IGA iga,IGAUserIFunction  IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIJacobian (IGA iga,IGAUserIJacobian  IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIFunction2(IGA iga,IGAUserIFunction2 IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIJacobian2(IGA iga,IGAUserIJacobian2 IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *ctx);
 
 PETSC_EXTERN PetscErrorCode IGAGetMeshInformation(IGA iga,PetscReal *hmin,PetscReal *hmax,PetscReal *havg,PetscReal *hstd);
 
   PetscInt  width[3];
   PetscInt  sizes[3];
   PetscInt  ID[3];
+  PetscBool atboundary;
+  PetscInt  atboundary_id;
   /**/
   PetscInt count;
   PetscInt index;
                        /*1: [nqp][nen][nsd]           */
                        /*2: [nqp][nen][nsd][nsd]      */
                        /*3: [nqp][nen][nsd][nsd][nsd] */
-  
+
   PetscReal *normal;   /*   [nqp][dim]                */
-  
+
 
   IGA      parent;
   IGAPoint iterator;
 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 PetscBool      IGAElementNextUserOps(IGAElement element,IGAUserOps *userops);
 PETSC_EXTERN PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAElementBeginPoint(IGAElement element,IGAPoint *point);
 PETSC_EXTERN PetscBool      IGAElementNextPoint(IGAElement element,IGAPoint point);
 PETSC_EXTERN PetscErrorCode IGAElementBuildProperty(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildQuadrature(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildShapeFuns(IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAElementBuildQuadratureAtBoundary(IGAElement element,PetscInt dir,PetscInt side);
+PETSC_EXTERN PetscErrorCode IGAElementBuildShapeFunsAtBoundary(IGAElement element,PetscInt dir,PetscInt side);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetParent(IGAElement element,IGA *parent);
 PETSC_EXTERN PetscErrorCode IGAElementGetIndex(IGAElement element,PetscInt *index);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *neq,PetscInt *nen,PetscInt *dof);
 PETSC_EXTERN PetscErrorCode IGAElementGetMapping(IGAElement element,PetscInt *nen,const PetscInt *mapping[]);
-PETSC_EXTERN PetscErrorCode IGAElementGetQuadrature(IGAElement element,PetscInt *nqp,PetscInt *dim,
-                                                    const PetscReal *point[],const PetscReal *weigth[],
-                                                    const PetscReal *detJac[]);
-PETSC_EXTERN PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,
-                                                   PetscInt *nen,PetscInt *dim,
-                                                   const PetscReal **shapefuns[]);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetWorkVal(IGAElement element,PetscScalar *U[]);
 PETSC_EXTERN PetscErrorCode IGAElementGetWorkVec(IGAElement element,PetscScalar *V[]);
 PETSC_EXTERN PetscErrorCode IGAElementAssembleVec(IGAElement element,const PetscScalar F[],Vec vec);
 PETSC_EXTERN PetscErrorCode IGAElementAssembleMat(IGAElement element,const PetscScalar K[],Mat mat);
 
-PETSC_EXTERN PetscBool      IGANextBoundaryElement(IGA iga,IGAElement element,PetscInt dir,PetscInt side);
-PETSC_EXTERN PetscErrorCode IGABoundaryElementBeginPoint(IGAElement element,IGAPoint *point,PetscInt dir,PetscInt side);
-PETSC_EXTERN PetscErrorCode IGAElementBuildBoundaryQuadrature(IGAElement element,PetscInt dir,PetscInt side);
-PETSC_EXTERN PetscErrorCode IGAElementBuildBoundaryShapeFuns(IGAElement element,PetscInt dir,PetscInt side);
-
 PETSC_EXTERN PetscErrorCode IGAElementCharacteristicSize(IGAElement element,PetscReal *h);
 
 /* ---------------------------------------------------------------- */
 
   PetscReal *weight;   /*   [1]   */
   PetscReal *detJac;   /*   [1]   */
+
   PetscReal *point;    /*   [dim] */
   PetscReal *scale;    /*   [dim] */
   PetscReal *basis[4]; /*0: [nen] */
 #define PCBBB "bbb"
 PETSC_EXTERN PetscErrorCode IGACreateKSP(IGA iga,KSP *ksp);
 PETSC_EXTERN PetscErrorCode IGAComputeSystem(IGA iga,Mat A,Vec B);
-PETSC_EXTERN PetscErrorCode IGAFormSystem(IGA iga,Mat A,Vec B,
-                                          IGAUserSystem,void *);
-PETSC_EXTERN PetscErrorCode IGABoundaryFormSystem(IGA iga,PetscInt dir,PetscInt side,
-						  Mat matA,Vec vecB,
-						  IGAUserSystem System,void *ctx);
 
 PETSC_EXTERN PetscErrorCode IGACreateSNES(IGA iga,SNES *snes);
-PETSC_EXTERN PetscErrorCode IGAFormFunction(IGA iga,Vec U,Vec F,IGAUserFunction,void *);
-PETSC_EXTERN PetscErrorCode IGAFormJacobian(IGA iga,Vec U,Mat J,IGAUserJacobian,void *);
 PETSC_EXTERN PetscErrorCode IGAComputeFunction(IGA iga,Vec U,Vec F);
 PETSC_EXTERN PetscErrorCode IGAComputeJacobian(IGA iga,Vec U,Mat J);
 
 PETSC_EXTERN PetscErrorCode IGACreateTS(IGA iga,TS *ts);
-PETSC_EXTERN PetscErrorCode IGAFormIFunction(IGA iga,PetscReal dt,
-                                             PetscReal a,Vec V,
-                                             PetscReal t,Vec U,
-                                             Vec F,IGAUserIFunction,void *);
-PETSC_EXTERN PetscErrorCode IGAFormIJacobian(IGA iga,PetscReal dt,
-                                             PetscReal a,Vec V,
-                                             PetscReal t,Vec U,
-                                             Mat J,IGAUserIJacobian,void *);
 PETSC_EXTERN PetscErrorCode IGAComputeIFunction(IGA iga,PetscReal dt,
                                                 PetscReal a,Vec V,
                                                 PetscReal t,Vec U,
                                                 PetscReal a,Vec V,
                                                 PetscReal t,Vec U,
                                                 Mat J);
-PETSC_EXTERN PetscErrorCode IGAFormIEFunction(IGA iga,PetscReal dt,
-                                              PetscReal a, Vec V,
-                                              PetscReal t, Vec U,
-                                              PetscReal t0,Vec U0,
-                                              Vec F,IGAUserIEFunction,void *);
-PETSC_EXTERN PetscErrorCode IGAFormIEJacobian(IGA iga,PetscReal dt,
-                                              PetscReal a, Vec V,
-                                              PetscReal t, Vec U,
-                                              PetscReal t0,Vec U0,
-                                              Mat J,IGAUserIEJacobian,void *);
 PETSC_EXTERN PetscErrorCode IGAComputeIEFunction(IGA iga,PetscReal dt,
                                                  PetscReal a, Vec V,
                                                  PetscReal t, Vec U,
                                                  Mat J);
 
 PETSC_EXTERN PetscErrorCode IGACreateTS2(IGA iga, TS *ts);
-PETSC_EXTERN PetscErrorCode IGAFormIFunction2(IGA iga,PetscReal dt,
-                                              PetscReal a,Vec vecA,
-                                              PetscReal v,Vec vecV,
-                                              PetscReal t,Vec vecU,
-                                              Vec vecF,
-                                              IGAUserIFunction2 IFunction, void *ctx);
-PETSC_EXTERN PetscErrorCode IGAFormIJacobian2(IGA iga,PetscReal dt,
-                                              PetscReal a,Vec vecA,
-                                              PetscReal v,Vec vecV,
-                                              PetscReal t,Vec vecU,
-                                              Mat matJ,
-                                              IGAUserIJacobian2 IJacobian,void *ctx);
 PETSC_EXTERN PetscErrorCode IGAComputeIFunction2(IGA iga,PetscReal dt,
                                                  PetscReal a,Vec vecA,
                                                  PetscReal v,Vec vecV,

src/petigabasis.c

   basis->value  = value;
 
   {
-    PetscInt  k0 = span[0], k1 = span[nel-1];
-    PetscReal u0 = U[k0],   u1 = U[k1+1];
+    PetscInt  o0 = offset[0], o1 = offset[nel-1];
+    PetscInt  k0 = span[0],   k1 = span[nel-1];
+    PetscReal u0 = U[k0],     u1 = U[k1+1];
     ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[0]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[1]);CHKERRQ(ierr);
+    basis->bnd_offset[0] =  o0; basis->bnd_offset[1] =  o1;
     basis->bnd_detJ  [0] = 1.0; basis->bnd_detJ  [1] = 1.0;
     basis->bnd_weight[0] = 1.0; basis->bnd_weight[1] = 1.0;
-    basis->bnd_point [0] =  u0; basis->bnd_point [1] = u1;
+    basis->bnd_point [0] =  u0; basis->bnd_point [1] =  u1;
     IGA_Basis_BSpline(k0,u0,p,d,U,basis->bnd_value[0]);
     IGA_Basis_BSpline(k1,u1,p,d,U,basis->bnd_value[1]);
   }

src/petigabound.c

 #define __FUNCT__ "IGABoundaryClear"
 PetscErrorCode IGABoundaryClear(IGABoundary boundary)
 {
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(boundary,1);
-  boundary->count   = 0;
+  boundary->count = 0;
   boundary->nload = 0;
+  ierr = PetscMemzero(boundary->userops,sizeof(struct _IGAUserOps));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   if (SysCtx) boundary->userops->SysCtx = SysCtx;
   PetscFunctionReturn(0);
 }
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserFunction"
+PetscErrorCode IGABoundarySetUserFunction(IGABoundary boundary,IGAUserFunction Function,void *FunCtx)
+{
+  PetscFunctionBegin;
+  if (Function) boundary->userops->Function = Function;
+  if (FunCtx)   boundary->userops->FunCtx   = FunCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserJacobian"
+PetscErrorCode IGABoundarySetUserJacobian(IGABoundary boundary,IGAUserJacobian Jacobian,void *JacCtx)
+{
+  PetscFunctionBegin;
+  if (Jacobian) boundary->userops->Jacobian = Jacobian;
+  if (JacCtx)   boundary->userops->JacCtx   = JacCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIFunction"
+PetscErrorCode IGABoundarySetUserIFunction(IGABoundary boundary,IGAUserIFunction IFunction,void *IFunCtx)
+{
+  PetscFunctionBegin;
+  if (IFunction) boundary->userops->IFunction = IFunction;
+  if (IFunCtx)   boundary->userops->IFunCtx   = IFunCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIJacobian"
+PetscErrorCode IGABoundarySetUserIJacobian(IGABoundary boundary,IGAUserIJacobian IJacobian,void *IJacCtx)
+{
+  PetscFunctionBegin;
+  if (IJacobian) boundary->userops->IJacobian = IJacobian;
+  if (IJacCtx)   boundary->userops->IJacCtx   = IJacCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIFunction2"
+PetscErrorCode IGABoundarySetUserIFunction2(IGABoundary boundary,IGAUserIFunction2 IFunction,void *IFunCtx)
+{
+  PetscFunctionBegin;
+  if (IFunction) boundary->userops->IFunction2 = IFunction;
+  if (IFunCtx)   boundary->userops->IFunCtx    = IFunCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIJacobian2"
+PetscErrorCode IGABoundarySetUserIJacobian2(IGABoundary boundary,IGAUserIJacobian2 IJacobian,void *IJacCtx)
+{
+  PetscFunctionBegin;
+  if (IJacobian) boundary->userops->IJacobian2 = IJacobian;
+  if (IJacCtx)   boundary->userops->IJacCtx    = IJacCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIEFunction"
+PetscErrorCode IGABoundarySetUserIEFunction(IGABoundary boundary,IGAUserIEFunction IEFunction,void *IEFunCtx)
+{
+  PetscFunctionBegin;
+  if (IEFunction) boundary->userops->IEFunction = IEFunction;
+  if (IEFunCtx)   boundary->userops->IEFunCtx   = IEFunCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserIEJacobian"
+PetscErrorCode IGABoundarySetUserIEJacobian(IGABoundary boundary,IGAUserIEJacobian IEJacobian,void *IEJacCtx)
+{
+  PetscFunctionBegin;
+  if (IEJacobian) boundary->userops->IEJacobian = IEJacobian;
+  if (IEJacCtx)   boundary->userops->IEJacCtx   = IEJacCtx;
+  PetscFunctionReturn(0);
+}
     ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&element->shape[2]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&element->shape[3]);CHKERRQ(ierr);
 
+    ierr = PetscMalloc1(nqp*dim,PetscReal,&element->normal);CHKERRQ(ierr);
+
+    /* */
 
     ierr = PetscMemzero(element->weight,  sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     ierr = PetscMemzero(element->detJac,  sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     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 = PetscMalloc1(nqp*dim,PetscReal,&element->normal);CHKERRQ(ierr);
     ierr = PetscMemzero(element->normal,sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
   }
   if (!element->collocation) {
     element->index = -1;
     return PETSC_FALSE;
   }
-  /* */
   for (i=0; i<dim; i++) {
     coord = index % width[i];
     index = (index - coord) / width[i];
     ID[i] = coord + start[i];
   }
-#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 = IGAElementBuildProperty(element); CHKERRRETURN(ierr,PETSC_FALSE);
-  ierr = IGAElementBuildFix(element);      CHKERRRETURN(ierr,PETSC_FALSE);
-#undef  CHKERRRETURN
-  return PETSC_TRUE;
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGANextBoundaryElement"
-PetscBool IGANextBoundaryElement(IGA iga,IGAElement element,PetscInt dir,PetscInt side)
-{
-  PetscInt i,dim  = element->dim;
-  PetscInt *start = element->start;
-  PetscInt *width = element->width;
-  PetscInt *ID    = element->ID;
-  PetscInt nel    = iga->axis[dir]->nel;
-  PetscInt index,coord;
-  PetscErrorCode ierr;
   /* */
-  element->nval = 0;
-  element->nvec = 0;
-  element->nmat = 0;
-  /* */
-
-  do {
-    index = ++element->index;
-    if (PetscUnlikely(index >= element->count)) {
-      element->index = -1;
-      return PETSC_FALSE;
-    }
-    /* */
-    for (i=0; i<dim; i++) {
-      coord = index % width[i];
-      index = (index - coord) / width[i];
-      ID[i] = coord + start[i];
-    }
-  }while( ID[dir] != side*(nel-1) );
-
 #undef  CHKERRRETURN
 #define CHKERRRETURN(n,r) do{if(PetscUnlikely(n)){CHKERRCONTINUE(n);return(r);}}while(0)
   ierr = IGAElementBuildMapping(element);  CHKERRRETURN(ierr,PETSC_FALSE);
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAElementNextUserOps"
+PetscBool IGAElementNextUserOps(IGAElement element,IGAUserOps *userops)
+{
+  static PetscInt counter = 0; /* XXX */
+  IGA iga = element->parent;
+  PetscInt dim = element->dim;
+  for (; PetscUnlikely(counter < 2*dim); counter++) {
+    PetscInt dir  = counter / 2;
+    PetscInt side = counter % 2;
+    PetscInt iel  = side ? element->sizes[dir]-1 : 0;
+    if (element->ID[dir] != iel) continue;
+    element->atboundary = PETSC_TRUE;
+    element->atboundary_id = counter;
+    *userops = iga->boundary[dir][side]->userops;
+    if (!(*userops)) continue;
+    counter++;
+    return PETSC_TRUE;
+  }
+  if (PetscLikely(counter == 2*dim)) {
+    element->atboundary = PETSC_FALSE;
+    element->atboundary_id = 0;
+    *userops = iga->userops;
+    counter++;
+    return PETSC_TRUE;
+  }
+  *userops = 0;
+  counter = 0;
+  return PETSC_FALSE;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetPoint"
 PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point)
 {
   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;
   (*point)->count = element->nqp;
-  (*point)->index = -1;
 
   (*point)->neq = element->neq;
   (*point)->nen = element->nen;
   (*point)->nsd = element->nsd;
   (*point)->npd = element->npd;
 
-  PetscFunctionReturn(0);
-}
-
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGABoundaryElementBeginPoint"
-PetscErrorCode IGABoundaryElementBeginPoint(IGAElement element,IGAPoint *point,PetscInt dir,PetscInt side)
-{
-  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 = IGAElementBuildBoundaryQuadrature(element,dir,side);CHKERRQ(ierr);
-  ierr = IGAElementBuildBoundaryShapeFuns(element,dir,side);CHKERRQ(ierr);
-  *point = element->iterator;
-
-  (*point)->count = element->nqp / element->parent->rule[dir]->nqp;
-  (*point)->index = -1;
-
-  (*point)->neq = element->neq;
-  (*point)->nen = element->nen;
-  (*point)->dof = element->dof;
-  (*point)->dim = element->dim;
-  (*point)->nsd = element->nsd;
-  (*point)->npd = element->npd;
-
+  if (PetscUnlikely(!element->atboundary)) {
+    ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
+    ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
+  } else {
+    PetscInt i = element->atboundary_id / 2;
+    PetscInt s = element->atboundary_id % 2;
+    (*point)->count = element->nqp / element->BD[i]->nqp;
+    ierr = IGAElementBuildQuadratureAtBoundary(element,i,s);CHKERRQ(ierr);
+    ierr = IGAElementBuildShapeFunsAtBoundary (element,i,s);CHKERRQ(ierr);
+    PetscFunctionReturn(0);
+  }
   PetscFunctionReturn(0);
 }
 
   PetscFunctionReturn(0);
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementGetQuadrature"
-PetscErrorCode IGAElementGetQuadrature(IGAElement element,PetscInt *nqp,PetscInt *dim,
-				       const PetscReal *point[],const PetscReal *weight[],
-				       const PetscReal *detJac[])
-{
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  if (nqp)    PetscValidIntPointer(nqp,2);
-  if (dim)    PetscValidIntPointer(dim,3);
-  if (point)  PetscValidPointer(point,4);
-  if (weight) PetscValidPointer(weight,5);
-  if (detJac) PetscValidPointer(detJac,6);
-  if (nqp)    *nqp    = element->nqp;
-  if (dim)    *dim    = element->dim;
-  if (point)  *point  = element->point;
-  if (weight) *weight = element->weight;
-  if (detJac) *detJac = element->detJac;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementGetShapeFuns"
-PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,
-				      PetscInt *nen,PetscInt *dim,
-				      const PetscReal **shapefuns[])
-{
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  if (nqp)       PetscValidIntPointer(nqp,2);
-  if (nen)       PetscValidIntPointer(nen,3);
-  if (dim)       PetscValidIntPointer(dim,4);
-  if (shapefuns) PetscValidPointer(shapefuns,5);
-  if (nqp)       *nqp       = element->nqp;
-  if (nen)       *nqp       = element->nqp;
-  if (dim)       *dim       = element->dim;
-  if (shapefuns) *shapefuns = (const PetscReal **)element->shape;
-  PetscFunctionReturn(0);
-}
-
 /* -------------------------------------------------------------------------- */
 
 #undef  __FUNCT__
     PetscInt a=0, *mapping = element->mapping;
     for (ka=0; ka<knen; ka++) {
       for (ja=0; ja<jnen; ja++) {
-	for (ia=0; ia<inen; ia++) {
-	  PetscInt iA = (ioffset + ia) - istart;
-	  PetscInt jA = (joffset + ja) - jstart;
-	  PetscInt kA = (koffset + ka) - kstart;
-	  mapping[a++] = iA + jA*jstride + kA*kstride;
-	}
+        for (ia=0; ia<inen; ia++) {
+          PetscInt iA = (ioffset + ia) - istart;
+          PetscInt jA = (joffset + ja) - jstart;
+          PetscInt kA = (koffset + ka) - kstart;
+          mapping[a++] = iA + jA*jstride + kA*kstride;
+        }
       }
     }
   }
     PetscInt i,nsd = element->nsd;
     if (element->geometry)
       for (a=0; a<nen; a++)
-	for (i=0; i<nsd; i++)
-	  X[i + a*nsd] = arrayX[i + map[a]*nsd];
+        for (i=0; i<nsd; i++)
+          X[i + a*nsd] = arrayX[i + map[a]*nsd];
     if (element->rational)
       for (a=0; a<nen; a++)
-	W[a] = arrayW[map[a]];
+        W[a] = arrayW[map[a]];
   }
   PetscFunctionReturn(0);
 }
     PetscInt i,npd = element->npd;
     for (a=0; a<nen; a++)
       for (i=0; i<npd; i++)
-	A[i + a*npd] = arrayA[i + map[a]*npd];
+        A[i + a*npd] = arrayA[i + map[a]*npd];
   }
   PetscFunctionReturn(0);
 }
 
+/* -------------------------------------------------------------------------- */
+
 EXTERN_C_BEGIN
 extern void IGA_Quadrature_1D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
+                              PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
 extern void IGA_Quadrature_2D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
+                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
+                              PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
 extern void IGA_Quadrature_3D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
-			      PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
+                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
+                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
+                              PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
+EXTERN_C_END
+
+EXTERN_C_BEGIN
+extern void IGA_BasisFuns_1D(PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+extern void IGA_BasisFuns_2D(PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+extern void IGA_BasisFuns_3D(PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+EXTERN_C_END
+
+EXTERN_C_BEGIN
+extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[]);
+extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[]);
+extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
+                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],
+                             PetscReal[],PetscReal[],PetscReal[]);
 EXTERN_C_END
 
 #define IGA_Quadrature_ARGS(ID,BD,i) \
   BD[i]->nqp,BD[i]->point+ID[i]*BD[i]->nqp,BD[i]->weight,BD[i]->detJ+ID[i]
 
+#define IGA_BasisFuns_ARGS(ID,BD,i) \
+  BD[i]->nqp,BD[i]->nen,BD[i]->d,BD[i]->value+ID[i]*BD[i]->nqp*BD[i]->nen*(BD[i]->d+1)
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementBuildQuadrature"
 PetscErrorCode IGAElementBuildQuadrature(IGAElement element)
   {
     IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
+    PetscReal *w = element->weight;
+    PetscReal *J = element->detJac;
+    PetscReal *u = element->point;
+    PetscReal *L = element->scale;
     switch (element->dim) {
     case 3: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
-			      IGA_Quadrature_ARGS(ID,BD,1),
-			      IGA_Quadrature_ARGS(ID,BD,2),
-			      element->weight,element->detJac,
-			      element->point,element->scale); break;
+                              IGA_Quadrature_ARGS(ID,BD,1),
+                              IGA_Quadrature_ARGS(ID,BD,2),
+                              w,J,u,L); break;
     case 2: IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
-			      IGA_Quadrature_ARGS(ID,BD,1),
-			      element->weight,element->detJac,
-			      element->point,element->scale); break;
+                              IGA_Quadrature_ARGS(ID,BD,1),
+                              w,J,u,L); break;
     case 1: IGA_Quadrature_1D(IGA_Quadrature_ARGS(ID,BD,0),
-			      element->weight,element->detJac,
-			      element->point,element->scale); break;
+                              w,J,u,L); break;
     }
   }
   PetscFunctionReturn(0);
 }
 
-#define IGA_Quadrature_BNDS(ID,BD,i,s) \
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementBuildShapeFuns"
+PetscErrorCode IGAElementBuildShapeFuns(IGAElement element)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  {
+    IGABasis  *BD = element->BD;
+    PetscInt  *ID = element->ID;
+    PetscInt  ord = element->parent->order;
+    PetscInt  rat = element->rational;
+    PetscReal *W  = element->rationalW;
+    PetscReal **N = element->basis;
+    switch (element->dim) {
+    case 3: IGA_BasisFuns_3D(ord,rat,W,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             IGA_BasisFuns_ARGS(ID,BD,1),
+                             IGA_BasisFuns_ARGS(ID,BD,2),
+                             N[0],N[1],N[2],N[3]); break;
+    case 2: IGA_BasisFuns_2D(ord,rat,W,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             IGA_BasisFuns_ARGS(ID,BD,1),
+                             N[0],N[1],N[2],N[3]); break;
+    case 1: IGA_BasisFuns_1D(ord,rat,W,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             N[0],N[1],N[2],N[3]); break;
+    }
+  }
+
+  if (element->dim == element->nsd) /* XXX */
+  if (element->geometry) {
+    PetscInt q;
+    PetscInt ord  = element->parent->order;
+    PetscInt nqp  = element->nqp;
+    PetscInt nen  = element->nen;
+    PetscReal *X  = element->geometryX;
+    PetscReal **M = element->basis;
+    PetscReal **N = element->shape;
+    PetscReal *J  = element->detX;
+    PetscReal *F  = element->gradX[0];
+    PetscReal *G  = element->gradX[1];
+    switch (element->dim) {
+    case 3: IGA_ShapeFuns_3D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
+    case 2: IGA_ShapeFuns_2D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
+    case 1: IGA_ShapeFuns_1D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
+    }
+    for (q=0; q<nqp; q++)
+      element->detJac[q] *= J[q];
+  }
+  PetscFunctionReturn(0);
+}
+
+#define IGA_Quadrature_BNDR(ID,BD,i,s) \
   1,&BD[i]->bnd_point[s],&BD[i]->bnd_weight[s],&BD[i]->bnd_detJ[s]
 
+#define IGA_BasisFuns_BNDR(ID,BD,i,s) \
+  1,BD[i]->nen,BD[i]->d,BD[i]->bnd_value[s]
+
+EXTERN_C_BEGIN
+extern void IGA_GetNormal(PetscInt dim,PetscInt dir,PetscInt side,const PetscReal F[],PetscReal *dS,PetscReal N[]);
+EXTERN_C_END
+
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementBuildBoundaryQuadrature"
-PetscErrorCode IGAElementBuildBoundaryQuadrature(IGAElement element,PetscInt dir,PetscInt side)
+#define __FUNCT__ "IGAElementBuildQuadratureAtBoundary"
+PetscErrorCode IGAElementBuildQuadratureAtBoundary(IGAElement element,PetscInt dir,PetscInt side)
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   {
     IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
+    PetscReal *w = element->weight;
+    PetscReal *J = element->detJac;
+    PetscReal *u = element->point;
+    PetscReal *L = element->scale;
     switch (element->dim) {
     case 3:
       switch (dir) {
-      case 0: IGA_Quadrature_3D(IGA_Quadrature_BNDS(ID,BD,0,side),
+      case 0: IGA_Quadrature_3D(IGA_Quadrature_BNDR(ID,BD,0,side),
                                 IGA_Quadrature_ARGS(ID,BD,1),
                                 IGA_Quadrature_ARGS(ID,BD,2),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
+                                w,J,u,L); break;
       case 1: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
-                                IGA_Quadrature_BNDS(ID,BD,1,side),
+                                IGA_Quadrature_BNDR(ID,BD,1,side),
                                 IGA_Quadrature_ARGS(ID,BD,2),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
+                                w,J,u,L); break;
       case 2: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
                                 IGA_Quadrature_ARGS(ID,BD,1),
-                                IGA_Quadrature_BNDS(ID,BD,2,side),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
-      }
-      break;
+                                IGA_Quadrature_BNDR(ID,BD,2,side),
+                                w,J,u,L); break;
+      } break;
     case 2:
       switch (dir) {
-      case 0: IGA_Quadrature_2D(IGA_Quadrature_BNDS(ID,BD,0,side),
+      case 0: IGA_Quadrature_2D(IGA_Quadrature_BNDR(ID,BD,0,side),
                                 IGA_Quadrature_ARGS(ID,BD,1),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
+                                w,J,u,L); break;
       case 1: IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
-                                IGA_Quadrature_BNDS(ID,BD,1,side),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
-      }
-      break;
+                                IGA_Quadrature_BNDR(ID,BD,1,side),
+                                w,J,u,L); break;
+      } break;
     case 1:
       switch (dir) {
-      case 0: IGA_Quadrature_1D(IGA_Quadrature_BNDS(ID,BD,0,side),
-                                element->weight,element->detJac,
-                                element->point,element->scale); break;
-      }
-      break;
+      case 0: IGA_Quadrature_1D(IGA_Quadrature_BNDR(ID,BD,0,side),
+                                w,J,u,L); break;
+      } break;
     }
   }
   PetscFunctionReturn(0);
 }
 
-EXTERN_C_BEGIN
-extern void IGA_BasisFuns_1D(PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
-extern void IGA_BasisFuns_2D(PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
-extern void IGA_BasisFuns_3D(PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
-EXTERN_C_END
-
-EXTERN_C_BEGIN
-extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
-extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
-extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
-			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
-EXTERN_C_END
-
-#define IGA_BasisFuns_ARGS(ID,BD,i) \
-  BD[i]->nqp,BD[i]->nen,BD[i]->d,   \
-  BD[i]->value+ID[i]*BD[i]->nqp*BD[i]->nen*(BD[i]->d+1)
-
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementBuildShapeFuns"
-PetscErrorCode IGAElementBuildShapeFuns(IGAElement element)
+#define __FUNCT__ "IGAElementBuildShapeFunsAtBoundary"
+PetscErrorCode IGAElementBuildShapeFunsAtBoundary(IGAElement element,PetscInt dir,PetscInt side)
 {
-  PetscInt order;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  order = element->parent->order;
   {
-    IGABasis *BD  = element->BD;
-    PetscInt *ID  = element->ID;
-    PetscReal **N = element->basis;
-    switch (element->dim) {
-    case 3: IGA_BasisFuns_3D(order,element->rational,
-			     element->rationalW,
-			     IGA_BasisFuns_ARGS(ID,BD,0),
-			     IGA_BasisFuns_ARGS(ID,BD,1),
-			     IGA_BasisFuns_ARGS(ID,BD,2),
-			     N[0],N[1],N[2],N[3]); break;
-    case 2: IGA_BasisFuns_2D(order,element->rational,
-			     element->rationalW,
-			     IGA_BasisFuns_ARGS(ID,BD,0),
-			     IGA_BasisFuns_ARGS(ID,BD,1),
-			     N[0],N[1],N[2],N[3]); break;
-    case 1: IGA_BasisFuns_1D(order,element->rational,
-			     element->rationalW,
-			     IGA_BasisFuns_ARGS(ID,BD,0),
-			     N[0],N[1],N[2],N[3]); break;
-    }
-  }
-  if (element->dim == element->nsd) /* XXX */
-  if (element->geometry) {
-    PetscInt q, nqp = element->nqp;
-    PetscReal **M = element->basis;
-    PetscReal **N = element->shape;
-    PetscReal *dX = element->detX;
-    PetscReal **gX = element->gradX;
-    switch (element->dim) {
-    case 3: IGA_ShapeFuns_3D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
-    case 2: IGA_ShapeFuns_2D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
-    case 1: IGA_ShapeFuns_1D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
-    }
-    for (q=0; q<nqp; q++)
-      element->detJac[q] *= dX[q];
-  }
-  PetscFunctionReturn(0);
-}
-
-#define IGA_BasisFuns_BNDS(ID,BD,i,s) \
-  1,BD[i]->nen,BD[i]->d,BD[i]->bnd_value[s]
-
-EXTERN_C_BEGIN
-extern void IGA_GetNormal(PetscInt dim,PetscInt dir,PetscInt side,const PetscReal F[],PetscReal *dS,PetscReal N[]);
-EXTERN_C_END
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementBuildBoundaryShapeFuns"
-PetscErrorCode IGAElementBuildBoundaryShapeFuns(IGAElement element,PetscInt dir,PetscInt side)
-{
-  PetscInt order;
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  if (PetscUnlikely(element->index < 0))
-    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  order = element->parent->order;
-  {
-    IGABasis *BD  = element->BD;
-    PetscInt *ID  = element->ID;
+    IGABasis  *BD = element->BD;
+    PetscInt  *ID = element->ID;
+    PetscInt  ord = element->parent->order;
+    PetscInt  rat = element->rational;
+    PetscReal *W  = element->rationalW;
     PetscReal **N = element->basis;
     switch (element->dim) {
     case 3:
       switch (dir) {
-      case 0: IGA_BasisFuns_3D(order,
-                               element->rational, element->rationalW,
-                               IGA_BasisFuns_BNDS(ID,BD,0,side),
+      case 0: IGA_BasisFuns_3D(ord,rat,W,
+                               IGA_BasisFuns_BNDR(ID,BD,0,side),
                                IGA_BasisFuns_ARGS(ID,BD,1),
                                IGA_BasisFuns_ARGS(ID,BD,2),
                                N[0],N[1],N[2],N[3]); break;
-      case 1: IGA_BasisFuns_3D(order,
-                               element->rational,element->rationalW,
+      case 1: IGA_BasisFuns_3D(ord,rat,W,
                                IGA_BasisFuns_ARGS(ID,BD,0),
-                               IGA_BasisFuns_BNDS(ID,BD,1,side),
+                               IGA_BasisFuns_BNDR(ID,BD,1,side),
                                IGA_BasisFuns_ARGS(ID,BD,2),
                                N[0],N[1],N[2],N[3]); break;
-      case 2: IGA_BasisFuns_3D(order,
-                               element->rational,element->rationalW,
+      case 2: IGA_BasisFuns_3D(ord,rat,W,
                                IGA_BasisFuns_ARGS(ID,BD,0),
                                IGA_BasisFuns_ARGS(ID,BD,1),
-                               IGA_BasisFuns_BNDS(ID,BD,2,side),
+                               IGA_BasisFuns_BNDR(ID,BD,2,side),
                                N[0],N[1],N[2],N[3]); break;
-      }
-      break;
+      } break;
     case 2:
       switch (dir) {
-      case 0: IGA_BasisFuns_2D(order,
-                               element->rational,element->rationalW,
-                               IGA_BasisFuns_BNDS(ID,BD,0,side),
+      case 0: IGA_BasisFuns_2D(ord,rat,W,
+                               IGA_BasisFuns_BNDR(ID,BD,0,side),
                                IGA_BasisFuns_ARGS(ID,BD,1),
                                N[0],N[1],N[2],N[3]); break;
-      case 1: IGA_BasisFuns_2D(order,
-                               element->rational,element->rationalW,
+      case 1: IGA_BasisFuns_2D(ord,rat,W,
                                IGA_BasisFuns_ARGS(ID,BD,0),
-                               IGA_BasisFuns_BNDS(ID,BD,1,side),
+                               IGA_BasisFuns_BNDR(ID,BD,1,side),
                                N[0],N[1],N[2],N[3]); break;
-      }
-      break;
-    case 1: 
+      } break;
+    case 1:
       switch (dir) {
-      case 0: IGA_BasisFuns_1D(order,
-                               element->rational,element->rationalW,
-                               IGA_BasisFuns_BNDS(ID,BD,0,side),
+      case 0: IGA_BasisFuns_1D(ord,rat,W,
+                               IGA_BasisFuns_BNDR(ID,BD,0,side),
                                N[0],N[1],N[2],N[3]); break;
-      }
-      break;
+      } break;
+    }
+    {
+      PetscInt q;
+      PetscInt nqp = element->nqp;
+      PetscInt dim = element->dim;
+      PetscReal *n = element->normal;
+      PetscMemzero(n,nqp*dim*sizeof(PetscReal));
+      for (q=0; q<nqp; q++)
+        n[q*dim+dir] = side ? 1.0 : -1.0;
     }
   }
+
   if (element->dim == element->nsd) /* XXX */
   if (element->geometry) {
-    PetscInt q, nqp = element->nqp;
+    PetscInt q;
+    PetscInt dim  = element->dim;
+    PetscInt ord  = element->parent->order;
+    PetscInt nqp  = element->nqp;
+    PetscInt nen  = element->nen;
+    PetscReal *X  = element->geometryX;
     PetscReal **M = element->basis;
     PetscReal **N = element->shape;
-    PetscReal *dX = element->detX;
-    PetscReal **gX = element->gradX;
-    switch (element->dim) {
-    case 3: IGA_ShapeFuns_3D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
-    case 2: IGA_ShapeFuns_2D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
-    case 1: IGA_ShapeFuns_1D(order,
-			     element->nqp,element->nen,
-			     element->geometryX,
-			     M[0],M[1],M[2],M[3],
-			     N[0],N[1],N[2],N[3],
-			     dX,gX[0],gX[1]); break;
+    PetscReal *J  = element->detX;
+    PetscReal *F  = element->gradX[0];
+    PetscReal *G  = element->gradX[1];
+    PetscReal *n  = element->normal;
+    switch (dim) {
+    case 3: IGA_ShapeFuns_3D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
+    case 2: IGA_ShapeFuns_2D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
+    case 1: IGA_ShapeFuns_1D(ord,nqp,nen,X,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             J,F,G); break;
     }
-    {
-      PetscInt dim = element->dim;
-      PetscReal *F = element->gradX[0];
-      PetscReal *n = element->normal;
-      for (q=0; q<nqp; q++) {
-        PetscReal dS;
-        IGA_GetNormal(dim,dir,side,&F[q*dim*dim],&dS,&n[q*dim]);
-        element->detJac[q] *= dS;
-      }
+    for (q=0; q<nqp; q++) {
+      PetscReal dS;
+      IGA_GetNormal(dim,dir,side,&F[q*dim*dim],&dS,&n[q*dim]);
+      element->detJac[q] *= dS;
     }
   }
-  if (!element->geometry) {
-    PetscInt q, nqp = element->nqp;
-    PetscInt i, dim = element->dim;
-    PetscReal *n = element->normal;
-    for (q=0; q<nqp; q++)
-      for (i=0; i<dim; i++)
-        n[q*dim+i] = (dir==i) ? ((side==0) ? -1.0 : 1.0) : 0.0;
-  }
   PetscFunctionReturn(0);
 }
 
     for (a=0; a<nen; a++) {
       const PetscScalar *u = arrayU + mapping[a]*dof;
       for (i=0; i<dof; i++)
-	U[pos++] = u[i]; /* XXX Use PetscMemcpy() ?? */
+        U[pos++] = u[i]; /* XXX Use PetscMemcpy() ?? */
     }
   }
   PetscFunctionReturn(0);
       PetscInt idx = a*dof + b->field[k];
       PetscScalar val = b->value[k];
       for (j=0; j<count; j++)
-	if (index[j] == idx) break;
+        if (index[j] == idx) break;
       if (j == count) count++;
       index[j] = idx;
       value[j] = val;
 
 EXTERN_C_BEGIN
 extern void IGA_BoundaryArea_2D(const PetscInt[],PetscInt,PetscInt,
-				PetscInt,const PetscReal[],
-				PetscInt,const PetscReal[],
-				PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
-				PetscReal*);
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscReal*);
 extern void IGA_BoundaryArea_3D(const PetscInt[],PetscInt,PetscInt,
-				PetscInt,const PetscReal[],
-				PetscInt,const PetscReal[],
-				PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
-				PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
-				PetscReal*);
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscReal*);
 EXTERN_C_END
 
 static PetscReal BoundaryArea(IGAElement element,PetscInt dir,PetscInt side)
     }
     switch (dim) {
     case 2: IGA_BoundaryArea_2D(shape,dir,side,
-				element->geometry,element->geometryX,
-				element->rational,element->rationalW,
-				nqp[0],W[0],nen[0],ndr[0],N[0],
-				&dS); break;
+                                element->geometry,element->geometryX,
+                                element->rational,element->rationalW,
+                                nqp[0],W[0],nen[0],ndr[0],N[0],
+                                &dS); break;
     case 3: IGA_BoundaryArea_3D(shape,dir,side,
-				element->geometry,element->geometryX,
-				element->rational,element->rationalW,
-				nqp[0],W[0],nen[0],ndr[0],N[0],
-				nqp[1],W[1],nen[1],ndr[1],N[1],
-				&dS);break;
+                                element->geometry,element->geometryX,
+                                element->rational,element->rationalW,
+                                nqp[0],W[0],nen[0],ndr[0],N[0],
+                                nqp[1],W[1],nen[1],ndr[1],N[1],
+                                &dS);break;
     }
     A *= dS;
   }
       PetscInt idx = a*dof + b->iload[k];
       PetscScalar val = b->vload[k];
       for (j=0; j<count; j++)
-	if (index[j] == idx) break;
+        if (index[j] == idx) break;
       if (j == count) value[count++] = 0.0;
       index[j] = idx;
       value[j]+= val*A;
     else      E[dir] = S[dir]+1;
     for (ka=S[2]; ka<E[2]; ka++)
       for (ja=S[1]; ja<E[1]; ja++)
-	for (ia=S[0]; ia<E[0]; ia++)
-	  {
-	    a = ia + ja*jstride + ka*kstride;
-	    AddFixa(element,b,a);
-	    AddFlux(element,b,a,Area);
-	  }
+        for (ia=S[0]; ia<E[0]; ia++)
+          {
+            a = ia + ja*jstride + ka*kstride;
+            AddFixa(element,b,a);
+            AddFlux(element,b,a,Area);
+          }
   }
 }
 
       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]) AddFixa(element,b[0][0],a);
-	      else if (iA == A1[0]) AddFixa(element,b[0][1],a);
-	      /**/ if (jA == A0[1]) AddFixa(element,b[1][0],a);
-	      else if (jA == A1[1]) AddFixa(element,b[1][1],a);
-	      /**/ if (kA == A0[2]) AddFixa(element,b[2][0],a);
-	      else if (kA == A1[2]) AddFixa(element,b[2][1],a);
-	      a++;
-	    }
+        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]) AddFixa(element,b[0][0],a);
+              else if (iA == A1[0]) AddFixa(element,b[0][1],a);
+              /**/ if (jA == A0[1]) AddFixa(element,b[1][0],a);
+              else if (jA == A1[1]) AddFixa(element,b[1][1],a);
+              /**/ if (kA == A0[2]) AddFixa(element,b[2][0],a);
+              else if (kA == A1[2]) AddFixa(element,b[2][1],a);
+              a++;
+            }
     }
   }
   PetscFunctionReturn(0);
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
       if (element->index == element->colmap[a]) {
-	PetscInt    c = k % element->dof;
-	PetscScalar v = element->vfix[f];
-	PetscInt i,j;
-	for (i=0; i<M; i++) K[i*N+k] = 0.0;
-	for (j=0; j<N; j++) K[c*N+j] = 0.0;
-	K[c*N+k] = 1.0;
-	F[c] = v;
+        PetscInt    c = k % element->dof;
+        PetscScalar v = element->vfix[f];
+        PetscInt i,j;
+        for (i=0; i<M; i++) K[i*N+k] = 0.0;
+        for (j=0; j<N; j++) K[c*N+j] = 0.0;
+        K[c*N+k] = 1.0;
+        F[c] = v;
       }
     }
   }
   PetscFunctionReturn(0);
 }
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementFixFunction"
 PetscErrorCode IGAElementFixFunction(IGAElement element,PetscScalar F[])
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
       if (element->index == 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;
+        PetscScalar v = element->vfix[f];
+        PetscScalar u = element->ufix[f];
+        F[c] = u - v;
       }
     }
   }
   PetscFunctionReturn(0);
 }
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementFixJacobian"
 PetscErrorCode IGAElementFixJacobian(IGAElement element,PetscScalar J[])
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
       if (element->index == element->colmap[a]) {
-	PetscInt c = k % element->dof;
-	PetscInt i,j;
-	if (element->index != element->colmap[a]) continue;
-	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 c = k % element->dof;
+        PetscInt i,j;
+        if (element->index != element->colmap[a]) continue;
+        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;
       }
     }
   }
   PetscFunctionReturn(0);
 }
 
+/* -------------------------------------------------------------------------- */
+
 EXTERN_C_BEGIN
 extern void IGA_Interpolate(PetscInt nen,PetscInt dof,PetscInt dim,PetscInt der,
-			    const PetscReal N[],const PetscScalar U[],PetscScalar u[]);
+                            const PetscReal N[],const PetscScalar U[],PetscScalar u[]);
 EXTERN_C_END
 
 #define ELENGTH(a,b) sqrt((corners[a][0]-corners[b][0])*(corners[a][0]-corners[b][0])+(corners[a][1]-corners[b][1])*(corners[a][1]-corners[b][1])+(corners[a][2]-corners[b][2])*(corners[a][2]-corners[b][2]))
 
 PETSC_EXTERN PetscLogEvent IGA_FormSystem;
 
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserSystem(IGAElement element,IGAUserSystem *sys,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->System);
+  if (!ops) return PETSC_FALSE;
+  *sys = ops->System;
+  *ctx = ops->SysCtx;
+  return PETSC_TRUE;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAComputeSystem"
 /*@
-   IGAFormSystem - Form the matrix and vector which represents the
+   IGAComputeSystem - Form the matrix and vector which represents the
    discretized a(w,u) = L(w).
 
    Collective on IGA/Mat/Vec
    Notes:
    This routine is used to solve a steady, linear problem. It performs
    matrix/vector assembly standard in FEM. The user provides a routine
-   which evaluates the bilinear form at a point.
+   which evaluates the bilinear and linear forms at a point.
 
    Level: normal
 
 @*/
 PetscErrorCode IGAComputeSystem(IGA iga,Mat matA,Vec vecB)
 {
+  IGAElement     element;
+  IGAPoint       point;
   IGAUserSystem  System;
-  void           *SysCtx;
-  PetscInt       dir,side;
+  void           *ctx;
+  PetscScalar    *A,*B;
+  PetscScalar    *K,*F;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
   IGACheckSetUp(iga,1);
   IGACheckUserOp(iga,1,System);
+
   ierr = MatZeroEntries(matA);CHKERRQ(ierr);
   ierr = VecZeroEntries(vecB);CHKERRQ(ierr);
-  for(dir=0; dir < iga->dim; dir++){
-    for(side=0; side < 2; side++){
-      System = iga->boundary[dir][side]->userops->System;
-      if(System){
-	SysCtx = iga->boundary[dir][side]->userops->SysCtx;
-	ierr = IGABoundaryFormSystem(iga,dir,side,matA,vecB,System,SysCtx);CHKERRQ(ierr);
-      }
-    }
-  }
-  System = iga->userops->System;
-  SysCtx = iga->userops->SysCtx;
-  ierr = IGAFormSystem(iga,matA,vecB,System,SysCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
 
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormSystem"
-PetscErrorCode IGAFormSystem(IGA iga,Mat matA,Vec vecB,
-			     IGAUserSystem System,void *ctx)
-{
-  IGAElement     element;
-  IGAPoint       point;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
-  PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
-  IGACheckSetUp(iga,1);
+  ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
 
   /* Element loop */
-  ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
   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);
-    /* 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 = System(point,K,F,ctx);CHKERRQ(ierr);
-      ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
-      ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
+    /* UserSystem loop */
+    while (IGAElementNextUserSystem(element,&System,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
+        ierr = IGAPointGetWorkVec(point,&F);CHKERRQ(ierr);
+        ierr = System(point,K,F,ctx);CHKERRQ(ierr);
+        ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
+        ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
+      }
+      ierr = IGAElementEndPoint(element,&point);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 = IGAEndElement(iga,&element);CHKERRQ(ierr);
+  
   ierr = PetscLogEventEnd(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
 
   ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGABoundaryFormSystem"
-PetscErrorCode IGABoundaryFormSystem(IGA iga,PetscInt dir,PetscInt side,
-				     Mat matA,Vec vecB,
-				     IGAUserSystem System,void *ctx)
-{
-  IGAElement     element;
-  IGAPoint       point;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
-  PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
-  IGACheckSetUp(iga,1);
-
-  /* Element loop */
-  ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
-  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
-  while (IGANextBoundaryElement(iga,element,dir,side)) {
-    PetscScalar *A, *B;
-    ierr = IGAElementGetWorkMat(element,&A);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&B);CHKERRQ(ierr);
-    /* Quadrature loop */
-    ierr = IGABoundaryElementBeginPoint(element,&point,dir,side);CHKERRQ(ierr);
-    while (IGAElementNextPoint(element,point)) {
-      PetscScalar *K, *F;
-      ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
-      ierr = IGAPointGetWorkVec(point,&F);CHKERRQ(ierr);
-      ierr = System(point,K,F,ctx);CHKERRQ(ierr);
-      ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
-      ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
-    }
-    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
-    ierr = IGAElementAssembleMat(element,A,matA);CHKERRQ(ierr);
-    ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
-  }
-  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
-  ierr = PetscLogEventEnd(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
-
-  ierr = MatAssemblyBegin(matA,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
-  ierr = MatAssemblyEnd  (matA,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
-  ierr = VecAssemblyBegin(vecB);CHKERRQ(ierr);
-  ierr = VecAssemblyEnd  (vecB);CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
 #define __FUNCT__ "IGACreateKSP"
 /*@
    IGACreateKSP - Creates a KSP (linear solver) which uses the same
 PETSC_EXTERN PetscLogEvent IGA_FormFunction;
 PETSC_EXTERN PetscLogEvent IGA_FormJacobian;
 
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserFunction(IGAElement element,IGAUserFunction *fun,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->Function);
+  if (!ops) return PETSC_FALSE;
+  *fun = ops->Function;
+  *ctx = ops->FunCtx;
+  return PETSC_TRUE;
+}
+
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserJacobian(IGAElement element,IGAUserJacobian *jac,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->Jacobian);
+  if (!ops) return PETSC_FALSE;
+  *jac = ops->Jacobian;
+  *ctx = ops->JacCtx;
+  return PETSC_TRUE;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAComputeFunction"
 PetscErrorCode IGAComputeFunction(IGA iga,Vec vecU,Vec vecF)
 {
+  Vec               localU;
+  const PetscScalar *arrayU;
+  IGAElement        element;
+  IGAPoint          point;
   IGAUserFunction   Function;
-  void              *FunCtx;
+  void              *ctx;
+  PetscScalar       *U,*F,*R;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(vecF,VEC_CLASSID,3);
   IGACheckSetUp(iga,1);
   IGACheckUserOp(iga,1,Function);
-  Function = iga->userops->Function;
-  FunCtx   = iga->userops->FunCtx;
-  ierr = IGAFormFunction(iga,vecU,vecF,Function,FunCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormFunction"
-PetscErrorCode IGAFormFunction(IGA iga,Vec vecU,Vec vecF,
-                               IGAUserFunction Function,void *ctx)
-{
-  Vec               localU;
-  const PetscScalar *arrayU;
-  IGAElement        element;
-  IGAPoint          point;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,2);
-  PetscValidHeaderSpecific(vecF,VEC_CLASSID,3);
-  IGACheckSetUp(iga,1);
 
   /* Clear global vector F*/
   ierr = VecZeroEntries(vecF);CHKERRQ(ierr);
   /* Get local vector U and array */
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecU,vecF,0);CHKERRQ(ierr);
+
   /* Element loop */
-  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecU,vecF,0);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *U, *F;
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
-    /* Quadrature loop */
-    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);
+    /* UserFunction loop */
+    while (IGAElementNextUserFunction(element,&Function,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
   ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
+
   ierr = PetscLogEventEnd(IGA_FormFunction,iga,vecU,vecF,0);CHKERRQ(ierr);
 
   /* Restore local vector U and array */
 #define __FUNCT__ "IGAComputeJacobian"
 PetscErrorCode IGAComputeJacobian(IGA iga,Vec vecU,Mat matJ)
 {
+  Vec               localU;
+  const PetscScalar *arrayU;
+  IGAElement        element;
+  IGAPoint          point;
   IGAUserJacobian   Jacobian;
-  void              *JacCtx;
+  void              *ctx;
+  PetscScalar       *U,*J,*K;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(matJ,MAT_CLASSID,3);
   IGACheckSetUp(iga,1);
   IGACheckUserOp(iga,1,Jacobian);
-  Jacobian = iga->userops->Jacobian;
-  JacCtx   = iga->userops->JacCtx;
-  ierr = IGAFormJacobian(iga,vecU,matJ,Jacobian,JacCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormJacobian"
-PetscErrorCode IGAFormJacobian(IGA iga,Vec vecU,Mat matJ,
-                               IGAUserJacobian Jacobian,void *ctx)
-{
-  Vec               localU;
-  const PetscScalar *arrayU;
-  IGAElement        element;
-  IGAPoint          point;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,2);
-  PetscValidHeaderSpecific(matJ,MAT_CLASSID,3);
-  IGACheckSetUp(iga,1);
 
   /* Clear global matrix J */
   ierr = MatZeroEntries(matJ);CHKERRQ(ierr);
   /* Get local vector U and array */
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecU,matJ,0);CHKERRQ(ierr);
+
   /* Element Loop */
-  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecU,matJ,0);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *U, *J;
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
     ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
     ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
-    /* Quadrature loop */
-    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);
+    /* UserJacobian loop */
+    while (IGAElementNextUserJacobian(element,&Jacobian,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
   }
   ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
+
   ierr = PetscLogEventEnd(IGA_FormJacobian,iga,vecU,matJ,0);CHKERRQ(ierr);
 
   /* Restore local vector U and array */
 
   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
   ierr = SNESCreate(comm,snes);CHKERRQ(ierr);
+  ierr = PetscObjectCompose((PetscObject)*snes,"IGA",(PetscObject)iga);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&F);CHKERRQ(ierr);
   ierr = SNESSetFunction(*snes,F,IGASNESFormFunction,iga);CHKERRQ(ierr);
   ierr = SNESSetJacobian(*snes,J,J,IGASNESFormJacobian,iga);CHKERRQ(ierr);
   ierr = MatDestroy(&J);CHKERRQ(ierr);
 
-  ierr = PetscObjectCompose((PetscObject)*snes,"IGA",(PetscObject)iga);CHKERRQ(ierr);
-
   PetscFunctionReturn(0);
 }
 PETSC_EXTERN PetscLogEvent IGA_FormFunction;
 PETSC_EXTERN PetscLogEvent IGA_FormJacobian;
 
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserIFunction(IGAElement element,IGAUserIFunction *fun,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->IFunction);
+  if (!ops) return PETSC_FALSE;
+  *fun = ops->IFunction;
+  *ctx = ops->IFunCtx;
+  return PETSC_TRUE;
+}
+
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserIJacobian(IGAElement element,IGAUserIJacobian *jac,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->IJacobian);
+  if (!ops) return PETSC_FALSE;
+  *jac = ops->IJacobian;
+  *ctx = ops->IJacCtx;
+  return PETSC_TRUE;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAComputeIFunction"
 PetscErrorCode IGAComputeIFunction(IGA iga,PetscReal dt,
                                    PetscReal t,Vec vecU,
                                    Vec vecF)
 {
-  IGAUserIFunction  IFunction;
-  void              *IFunCtx;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
-  PetscValidHeaderSpecific(vecF,VEC_CLASSID,7);
-  IGACheckSetUp(iga,1);
-  IGACheckUserOp(iga,1,IFunction);
-  IFunction = iga->userops->IFunction;
-  IFunCtx   = iga->userops->IFunCtx;
-  ierr = IGAFormIFunction(iga,dt,a,vecV,t,vecU,vecF,IFunction,IFunCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormIFunction"
-PetscErrorCode IGAFormIFunction(IGA iga,PetscReal dt,
-                                PetscReal a,Vec vecV,
-                                PetscReal t,Vec vecU,
-                                Vec vecF,
-                                IGAUserIFunction IFunction, void *ctx)
-{
   Vec               localV;
   Vec               localU;
   const PetscScalar *arrayV;
   const PetscScalar *arrayU;
   IGAElement        element;
   IGAPoint          point;
+  IGAUserIFunction  IFunction;
+  void              *ctx;
+  PetscScalar       *V,*U,*F,*R;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
   PetscValidHeaderSpecific(vecF,VEC_CLASSID,7);
   IGACheckSetUp(iga,1);
+  IGACheckUserOp(iga,1,IFunction);
 
   /* Clear global vector F */
   ierr = VecZeroEntries(vecF);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecV,&localV,&arrayV);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
+
   /* Element loop */
-  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *V, *U, *F;
     ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);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 = 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);
+    /* UserIFunction loop */
+    while (IGAElementNextUserIFunction(element,&IFunction,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);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 */
                                    PetscReal t,Vec vecU,
                                    Mat matJ)
 {
-  IGAUserIJacobian  IJacobian;
-  void              *IJacCtx;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
-  PetscValidHeaderSpecific(matJ,MAT_CLASSID,7);
-  IGACheckSetUp(iga,1);
-  IGACheckUserOp(iga,1,IJacobian);
-  IJacobian = iga->userops->IJacobian;
-  IJacCtx   = iga->userops->IJacCtx;
-  ierr = IGAFormIJacobian(iga,dt,a,vecV,t,vecU,matJ,IJacobian,IJacCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormIJacobian"
-PetscErrorCode IGAFormIJacobian(IGA iga,PetscReal dt,
-                                PetscReal a,Vec vecV,
-                                PetscReal t,Vec vecU,
-                                Mat matJ,
-                                IGAUserIJacobian IJacobian,void *ctx)
-{
   Vec               localV;
   Vec               localU;
   const PetscScalar *arrayV;
   const PetscScalar *arrayU;
   IGAElement        element;
   IGAPoint          point;
+  IGAUserIJacobian  IJacobian;
+  void              *ctx;
+  PetscScalar       *V,*U,*J,*K;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
+  PetscValidHeaderSpecific(vecV,VEC_CLASSID,4);
   PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
   PetscValidHeaderSpecific(matJ,MAT_CLASSID,7);
   IGACheckSetUp(iga,1);
+  IGACheckUserOp(iga,1,IJacobian);
 
   /* Clear global matrix J*/
   ierr = MatZeroEntries(matJ);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecV,&localV,&arrayV);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
+
   /* Element Loop */
-  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *V, *U, *J;
     ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);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 = 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);
+    /* UserIJacobian loop */
+    while (IGAElementNextUserIJacobian(element,&IJacobian,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);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 */
   PetscFunctionReturn(0);
 }
 
+
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserIEFunction(IGAElement element,IGAUserIEFunction *fun,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->IEFunction);
+  if (!ops) return PETSC_FALSE;
+  *fun = ops->IEFunction;
+  *ctx = ops->IEFunCtx;
+  return PETSC_TRUE;
+}
+
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserIEJacobian(IGAElement element,IGAUserIEJacobian *jac,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->IEJacobian);
+  if (!ops) return PETSC_FALSE;
+  *jac = ops->IEJacobian;
+  *ctx = ops->IEJacCtx;
+  return PETSC_TRUE;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAComputeIEFunction"
 PetscErrorCode IGAComputeIEFunction(IGA iga,PetscReal dt,
-                                    PetscReal a,Vec vecV,
-                                    PetscReal t,Vec vecU,
+                                    PetscReal a, Vec vecV,
+                                    PetscReal t, Vec vecU,
                                     PetscReal t0,Vec vecU0,
                                     Vec vecF)
 {
+  Vec               localV;
+  Vec               localU;
+  Vec               localU0;
+  const PetscScalar *arrayV;
+  const PetscScalar *arrayU;
+  const PetscScalar *arrayU0;
+  IGAElement        element;
+  IGAPoint          point;
   IGAUserIEFunction IEFunction;
-  void              *IEFunCtx;
+  void              *ctx;
+  PetscScalar       *V,*U,*U0,*F,*R;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(vecF,VEC_CLASSID,9);
   IGACheckSetUp(iga,1);
   IGACheckUserOp(iga,1,IEFunction);
-  IEFunction = iga->userops->IEFunction;
-  IEFunCtx   = iga->userops->IEFunCtx;
-  ierr = IGAFormIEFunction(iga,dt,a,vecV,t,vecU,t0,vecU0,vecF,IEFunction,IEFunCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormIEFunction"
-PetscErrorCode IGAFormIEFunction(IGA iga,PetscReal dt,
-                                 PetscReal a,Vec vecV,
-                                 PetscReal t,Vec vecU,
-                                 PetscReal t0,Vec vecU0,
-                                 Vec vecF,
-                                 IGAUserIEFunction IEFunction,void *ctx)
-{
-  Vec               localV;
-  Vec               localU;
-  Vec               localU0;
-  const PetscScalar *arrayV;
-  const PetscScalar *arrayU;
-  const PetscScalar *arrayU0;
-  IGAElement        element;
-  IGAPoint          point;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,4);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
-  PetscValidHeaderSpecific(vecU0,VEC_CLASSID,8);
-  PetscValidHeaderSpecific(vecF,VEC_CLASSID,9);
-  IGACheckSetUp(iga,1);
 
   /* Clear global vector F */
   ierr = VecZeroEntries(vecF);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU0,&localU0,&arrayU0);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
+
   /* Element loop */
-  ierr = PetscLogEventBegin(IGA_FormFunction,iga,vecV,vecU,vecF);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *V, *U, *U0, *F;
     ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U0);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);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);  /* XXX */
-    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
-    /* Quadrature loop */
-    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);
+    /* UserIEFunction loop */
+    while (IGAElementNextUserIEFunction(element,&IEFunction,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);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 */
 #undef  __FUNCT__
 #define __FUNCT__ "IGAComputeIEJacobian"
 PetscErrorCode IGAComputeIEJacobian(IGA iga,PetscReal dt,
-                                    PetscReal a,Vec vecV,
-                                    PetscReal t,Vec vecU,
+                                    PetscReal a, Vec vecV,
+                                    PetscReal t, Vec vecU,
                                     PetscReal t0,Vec vecU0,
                                     Mat matJ)
 {
-  IGAUserIEJacobian  IEJacobian;
-  void              *IEJacCtx;
-  PetscErrorCode    ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
-  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
-  PetscValidHeaderSpecific(matJ,MAT_CLASSID,7);
-  IGACheckSetUp(iga,1);
-  IGACheckUserOp(iga,1,IEJacobian);
-  IEJacobian = iga->userops->IEJacobian;
-  IEJacCtx   = iga->userops->IEJacCtx;
-  ierr = IGAFormIEJacobian(iga,dt,a,vecV,t,vecU,t0,vecU0,matJ,IEJacobian,IEJacCtx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAFormIEJacobian"
-PetscErrorCode IGAFormIEJacobian(IGA iga,PetscReal dt,
-                                 PetscReal a,Vec vecV,
-                                 PetscReal t,Vec vecU,
-                                 PetscReal t0,Vec vecU0,
-                                 Mat matJ,
-                                 IGAUserIEJacobian IEJacobian,void *ctx)
-{
   Vec               localV;
   Vec               localU;
   Vec               localU0;
   const PetscScalar *arrayU0;
   IGAElement        element;
   IGAPoint          point;
+  IGAUserIEJacobian IEJacobian;
+  void              *ctx;
+  PetscScalar       *V,*U,*U0,*J,*K;
   PetscErrorCode    ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
+  PetscValidHeaderSpecific(vecV,VEC_CLASSID,4);
   PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
-  PetscValidHeaderSpecific(vecU0,VEC_CLASSID,7);
-  PetscValidHeaderSpecific(matJ,MAT_CLASSID,0);
+  PetscValidHeaderSpecific(vecU0,VEC_CLASSID,8);
+  PetscValidHeaderSpecific(matJ,MAT_CLASSID,9);
   IGACheckSetUp(iga,1);
+  IGACheckUserOp(iga,1,IEJacobian);
 
   /* Clear global matrix J*/
   ierr = MatZeroEntries(matJ);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU,&localU,&arrayU);CHKERRQ(ierr);
   ierr = IGAGetLocalVecArray(iga,vecU0,&localU0,&arrayU0);CHKERRQ(ierr);
 
+  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
+
   /* Element Loop */
-  ierr = PetscLogEventBegin(IGA_FormJacobian,iga,vecV,vecU,matJ);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
   while (IGANextElement(iga,element)) {
-    PetscScalar *V, *U, *U0, *J;
     ierr = IGAElementGetWorkVal(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(element,&U0);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);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);  /* XXX */
-    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
-    /* Quadrature loop */
-    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);
+    /* UserIEJacobian loop */
+    while (IGAElementNextUserIEJacobian(element,&IEJacobian,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        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 = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
-    /* */
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);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 */
 {
   IGA            iga = (IGA)ctx;
   PetscReal      dt,a=0;
+  PetscReal      t0;
+  Vec            U0;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,6);
   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
   if (iga->userops->IEFunction) {
-    PetscReal t0;
-    Vec       U0;
     ierr = TSGetTime(ts,&t0);CHKERRQ(ierr);
     ierr = TSGetSolution(ts,&U0);CHKERRQ(ierr);
     ierr = IGAComputeIEFunction(iga,dt,a,V,t,U,t0,U0,F);CHKERRQ(ierr);
 {
   IGA            iga = (IGA)ctx;
   PetscReal      dt,a=shift;
+  PetscReal      t0;
+  Vec            U0;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,9);
   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
   if (iga->userops->IEJacobian) {
-    PetscReal t0;
-    Vec       U0;
     ierr = TSGetTime(ts,&t0);CHKERRQ(ierr);
     ierr = TSGetSolution(ts,&U0);CHKERRQ(ierr);
     ierr = IGAComputeIEJacobian(iga,dt,a,V,t,U,t0,U0,*P);CHKERRQ(ierr);
 
   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
   ierr = TSCreate(comm,ts);CHKERRQ(ierr);
+  ierr = PetscObjectCompose((PetscObject)*ts,"IGA",(PetscObject)iga);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = TSSetSolution(*ts,U);CHKERRQ(ierr);
   ierr = TSSetIJacobian(*ts,J,J,IGATSFormIJacobian,iga);CHKERRQ(ierr);
   ierr = MatDestroy(&J);CHKERRQ(ierr);
 
-  ierr = PetscObjectCompose((PetscObject)*ts,"IGA",(PetscObject)iga);CHKERRQ(ierr);