Commits

Nathan Collier committed f5201a1

added boundary integrals for linear problems

Comments (1)

  1. Nathan Collier author

    Take a look at how I implemented generic boundary integrals in our code framework. I only implemented the KSP part for the moment until we agree on how it should work. With this approach we can handle weak imposition of BCs, multipatch using DG-like fixes to merge, as well as general Neumann conditions. Also opens the door for fluid-structure interaction. I apologize for some wonky hacks--let me know what you think.

Files changed (6)

       PetscReal Nb_y = N1[b][1];
       K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
     }
+    F[a] = Na * 0.0;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Neumann"
+PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscReal *N0 = p->shape[0];
+  PetscInt a,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na   = N0[a];
     F[a] = Na * 1.0;
   }
   return 0;
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt N[2] = {16,16}, nN = 2; 
-  PetscInt p[2] = { 2, 2}, np = 2;
-  PetscInt C[2] = {-1,-1}, nC = 2;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&nN,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&np,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  if (nN == 1) N[1] = N[0];
-  if (np == 1) p[1] = p[0];
-  if (nC == 1) C[1] = C[0];
-  if (C[0] == -1) C[0] = p[0]-1;
-  if (C[1] == -1) C[1] = p[1]-1;
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
 
-  PetscInt i;
-  for (i=0; i<2; i++) {
-    IGAAxis axis;
-    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
-    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,N[i],-1.0,+1.0,C[i]);CHKERRQ(ierr);
-  }
   IGABoundary bnd;
-  PetscInt dir,side;
-  PetscScalar value = 1.0;
-  for (dir=0; dir<2; dir++) {
-    for (side=0; side<2; side++) {
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
-    }
-  }
-
+  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr); // u = 0 on [0,:]
+  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
+  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr); // grad u . n = h on [1,:]
+  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
 	  NavierStokesVMS \
 	  Elasticity3D \
 	  HyperElasticity \
-	  PhaseFieldCrystal2D
+	  PhaseFieldCrystal2D \
+          BoundaryIntegral
 
 ALL: ${TARGETS}
 clean::
 PhaseFieldCrystal2D: PhaseFieldCrystal2D.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
+BoundaryIntegral: BoundaryIntegral.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 
 OPTS=-nox -malloc_debug -malloc_dump
 
 PETSC_EXTERN PetscErrorCode IGABasisInitQuadrature (IGABasis basis,IGAAxis axis,IGARule rule,PetscInt order);
 PETSC_EXTERN PetscErrorCode IGABasisInitCollocation(IGABasis basis,IGAAxis axis,PetscInt order);
 
-struct _n_IGABoundary {
-  PetscInt refct;
-  /**/
-  PetscInt    dof;
-  /**/
-  PetscInt    count;
-  PetscInt    *field;
-  PetscScalar *value;
-  /**/
-  PetscInt    nload;
-  PetscInt    *iload;
-  PetscScalar *vload;
-};
-PETSC_EXTERN PetscErrorCode IGABoundaryCreate(IGABoundary *boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryDestroy(IGABoundary *boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryReset(IGABoundary boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryReference(IGABoundary boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryInit(IGABoundary boundary,PetscInt dof);
-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);
-
 /* ---------------------------------------------------------------- */
 
 typedef PetscErrorCode (*IGAUserScalar)    (IGAPoint point,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx);
 
 /* ---------------------------------------------------------------- */
 
+struct _n_IGABoundary {
+  PetscInt refct;
+  /**/
+  PetscInt    dof;
+  /**/
+  PetscInt    count;
+  PetscInt    *field;
+  PetscScalar *value;
+  /**/
+  PetscInt    nload;
+  PetscInt    *iload;
+  PetscScalar *vload;
+  /**/
+  IGAUserOps userops;
+};
+PETSC_EXTERN PetscErrorCode IGABoundaryCreate(IGABoundary *boundary);
+PETSC_EXTERN PetscErrorCode IGABoundaryDestroy(IGABoundary *boundary);
+PETSC_EXTERN PetscErrorCode IGABoundaryReset(IGABoundary boundary);
+PETSC_EXTERN PetscErrorCode IGABoundaryReference(IGABoundary boundary);
+PETSC_EXTERN PetscErrorCode IGABoundaryInit(IGABoundary boundary,PetscInt dof);
+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);
+
 typedef struct _p_IGA *IGA;
 
 typedef struct _IGAOps *IGAOps;
 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 PetscBool IGABoundaryElementNextPoint(IGAElement element,IGAPoint point,PetscInt dir,PetscInt side);
+
 /* ---------------------------------------------------------------- */
 
 struct _n_IGAPoint {
 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 *);

src/petigabound.c

   PetscFunctionBegin;
   PetscValidPointer(boundary,1);
   ierr = PetscNew(struct _n_IGABoundary,boundary);CHKERRQ(ierr);
+  ierr = PetscNew(struct _IGAUserOps,&((*boundary)->userops));CHKERRQ(ierr);
   (*boundary)->refct = 1;
   PetscFunctionReturn(0);
 }
   if (!boundary) PetscFunctionReturn(0);
   if (--boundary->refct > 0) PetscFunctionReturn(0);
   ierr = IGABoundaryReset(boundary);CHKERRQ(ierr);
+  ierr = PetscFree(boundary->userops);CHKERRQ(ierr);
   ierr = PetscFree(boundary);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /*@
    IGABoundarySetValue - Used to set a constant Dirichlet boundary
    condition on the given boundary.
-   
+
    Logically Collective on IGABoundary
 
    Input Parameters:
 /*@
    IGABoundarySetLoad - Used to set a constant Neumann boundary
    condition on the given boundary.
-   
+
    Logically Collective on IGABoundary
 
    Input Parameters:
   }
   PetscFunctionReturn(0);
 }
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserSystem"
+/*@
+   IGABoundarySetUserSystem - Set the user callback to form the matrix
+   and vector which represents the discretized a(w,u) = L(w)
+   integrated along the given boundary.
+
+   Logically collective on IGABoundary
+
+   Input Parameters:
++  boundary - the IGABoundary context
+.  System - the function which evaluates a(w,u) and L(w)
+-  ctx - user-defined context for evaluation routine (may be PETSC_NULL)
+
+   Details of System:
+$  PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx);
+
++  p - point at which to evaluate a(w,u)=L(w)
+.  K - contribution to a(w,u)
+.  F - contribution to L(w)
+-  ctx - user-defined context for evaluation routine
+
+   Level: normal
+
+.keywords: IGABoundary, setup linear system, matrix assembly, vector assembly
+@*/
+PetscErrorCode IGABoundarySetUserSystem(IGABoundary boundary,IGAUserSystem System,void *SysCtx)
+{
+  PetscFunctionBegin;
+  if (System) boundary->userops->System = System;
+  if (SysCtx) boundary->userops->SysCtx = SysCtx;
+  PetscFunctionReturn(0);
+}
 }
 
 #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);
+  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__ "IGAEndElement"
 PetscErrorCode IGAEndElement(IGA iga,IGAElement *element)
 {
   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;
+  (*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;
+
+  PetscFunctionReturn(0);
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementNextPoint"
 PetscBool IGAElementNextPoint(IGAElement element,IGAPoint point)
     point->geometry = PETSC_NULL;
   if (!element->property)
     point->property = PETSC_NULL;
-  
+
+  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__ "IGABoundaryElementNextPoint"
+PetscBool IGABoundaryElementNextPoint(IGAElement element,IGAPoint point,PetscInt dir,PetscInt side)
+{
+  PetscInt nen = point->nen;
+  PetscInt dim = point->dim;
+  PetscInt nsd = point->nsd;
+  PetscInt index,skip=1,i;
+  /* */
+  point->nvec = 0;
+  point->nmat = 0;
+  /* */
+  for(i=0;i<dim;i++) if(i != dir) skip *= element->parent->rule[i]->nqp;
+
+  do {
+    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;
+  } while( index % skip != 0 );
  1. Nathan Collier author

    This is wrong. It only works if you integrate along a x-direction boundary. Basically you went to great lengths to not loop over each dimension so we are dimension independent but I need something like it here.

+
+  return PETSC_TRUE;
+
+ start:
+
+  point->geometry = element->geometryX;
+  point->property = element->propertyA;
+  if (!element->geometry)
+    point->geometry = PETSC_NULL;
+  if (!element->property)
+    point->property = PETSC_NULL;
+
   point->weight   = element->weight;
   point->detJac   = element->detJac;
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetQuadrature"
 PetscErrorCode IGAElementGetQuadrature(IGAElement element,PetscInt *nqp,PetscInt *dim,
-                                       const PetscReal *point[],const PetscReal *weight[],
-                                       const PetscReal *detJac[])
+				       const PetscReal *point[],const PetscReal *weight[],
+				       const PetscReal *detJac[])
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetShapeFuns"
 PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,
-                                      PetscInt *nen,PetscInt *dim,
-                                      const PetscReal **shapefuns[])
+				      PetscInt *nen,PetscInt *dim,
+				      const PetscReal **shapefuns[])
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
     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
 
 #define IGA_Quadrature_ARGS(ID,BD,i) \
     PetscInt *ID = element->ID;
     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),
+			      element->weight,element->detJac,
+			      element->point,element->scale); 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),
+			      element->weight,element->detJac,
+			      element->point,element->scale); break;
     case 1: IGA_Quadrature_1D(IGA_Quadrature_ARGS(ID,BD,0),
-                              element->weight,element->detJac,
-                              element->point,element->scale); break;
+			      element->weight,element->detJac,
+			      element->point,element->scale); break;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementBuildBoundaryQuadrature"
+PetscErrorCode IGAElementBuildBoundaryQuadrature(IGAElement element,PetscInt dir,PetscInt side)
+{
+  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 i,*ID = element->ID;
+
+    /* Sad hack */
+    PetscScalar pts[BD[dir]->nqp],wts[BD[dir]->nqp];
+    PetscReal U[2] = {0,0}, one=1;
+    IGAAxisGetLimits(element->parent->axis[dir],&U[0],&U[1]);
+    for(i=0;i<BD[dir]->nqp;i++){
+      pts[i] = U[side];
+      wts[i] = 1.0;
+    }
+
+    switch (element->dim) {
+    case 3:
+      if(dir==0){
+	IGA_Quadrature_3D(BD[0]->nqp,&pts[0],&wts[0],&one,
+			  IGA_Quadrature_ARGS(ID,BD,1),
+			  IGA_Quadrature_ARGS(ID,BD,2),
+			  element->weight,element->detJac,
+			  element->point,element->scale);
+      }else if(dir==1){
+	IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
+			  BD[1]->nqp,&pts[0],&wts[0],&one,
+			  IGA_Quadrature_ARGS(ID,BD,2),
+			  element->weight,element->detJac,
+			  element->point,element->scale);
+      }else{
+	IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
+			  IGA_Quadrature_ARGS(ID,BD,1),
+			  BD[2]->nqp,&pts[0],&wts[0],&one,
+			  element->weight,element->detJac,
+			  element->point,element->scale);
+      }
+      break;
+    case 2:
+      if(dir==0){
+	IGA_Quadrature_2D(BD[0]->nqp,&pts[0],&wts[0],&one,
+			  IGA_Quadrature_ARGS(ID,BD,1),
+			  element->weight,element->detJac,
+			  element->point,element->scale);
+      }else{
+	IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
+			  BD[1]->nqp,&pts[0],&wts[0],&one,
+			  element->weight,element->detJac,
+			  element->point,element->scale);
+      }
+      break;
+    case 1: IGA_Quadrature_1D(IGA_Quadrature_ARGS(ID,BD,0),
+			      element->weight,element->detJac,
+			      element->point,element->scale); break;
     }
   }
   PetscFunctionReturn(0);
 
 EXTERN_C_BEGIN
 extern void IGA_BasisFuns_1D(PetscInt,PetscInt,const PetscReal[],
-                             PetscInt,PetscInt,PetscInt,const PetscReal[],
-                             PetscReal[],PetscReal[],PetscReal[],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[]);
+			     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[]);
+			     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[]);
+			     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[]);
+			     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[]);
+			     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) \
     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;
+			     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;
+			     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;
+			     element->rationalW,
+			     IGA_BasisFuns_ARGS(ID,BD,0),
+			     N[0],N[1],N[2],N[3]); break;
     }
   }
   if (element->dim == element->nsd) /* XXX */
     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;
+			     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;
+			     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;
+			     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);
+}
+
+extern void IGA_DersBasisFuns(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
+
+#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;
+    PetscReal **N = element->basis;
+
+    PetscInt k = element->parent->axis[dir]->span[element->ID[dir]];
+    PetscReal u = element->point[dir];
+    PetscInt p = element->parent->axis[dir]->p;
+    PetscInt d = element->parent->order;
+    PetscReal *U = element->parent->axis[dir]->U;
+    PetscInt nqp = BD[dir]->nqp;
+    PetscInt nen = BD[dir]->nen;
+    PetscReal value[nqp*nen*(d+1)];
+    IGA_DersBasisFuns(k,u,p,d,U,&value[0]);
+
+    switch (element->dim) {
+    case 3:
+      if(dir == 0){
+	IGA_BasisFuns_3D(order,element->rational,
+			 element->rationalW,
+			 BD[0]->nqp,BD[0]->nen,BD[0]->d,&value[0],
+			 IGA_BasisFuns_ARGS(ID,BD,1),
+			 IGA_BasisFuns_ARGS(ID,BD,2),
+			 N[0],N[1],N[2],N[3]);
+      }else if(dir == 1){
+	IGA_BasisFuns_3D(order,element->rational,
+			 element->rationalW,
+			 IGA_BasisFuns_ARGS(ID,BD,0),
+			 BD[1]->nqp,BD[1]->nen,BD[1]->d,&value[0],
+			 IGA_BasisFuns_ARGS(ID,BD,2),
+			 N[0],N[1],N[2],N[3]);
+      }else{
+	IGA_BasisFuns_3D(order,element->rational,
+			 element->rationalW,
+			 IGA_BasisFuns_ARGS(ID,BD,0),
+			 IGA_BasisFuns_ARGS(ID,BD,1),
+			 BD[2]->nqp,BD[2]->nen,BD[2]->d,&value[0],
+			 N[0],N[1],N[2],N[3]);
+      }
+      break;
+    case 2:
+      if(dir == 0){
+	IGA_BasisFuns_2D(order,element->rational,
+			 element->rationalW,
+			 BD[0]->nqp,BD[0]->nen,BD[0]->d,&value[0],
+			 IGA_BasisFuns_ARGS(ID,BD,1),
+			 N[0],N[1],N[2],N[3]);
+      }else{
+	IGA_BasisFuns_2D(order,element->rational,
+			 element->rationalW,
+			 IGA_BasisFuns_ARGS(ID,BD,0),
+			 BD[1]->nqp,BD[1]->nen,BD[1]->d,&value[0],
+			 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];
     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);
+	  }
   }
 }
 
       IGAAxis  *AX = element->parent->axis;
       PetscInt i,dim = element->dim;
       for (i=0; i<dim; i++) {
-        PetscBool w = AX[i]->periodic;
-        PetscInt  n = AX[i]->nnp-1; /* last node */
-        A0[i] = 0; if (!w) A1[i] = n;
+	PetscBool w = AX[i]->periodic;
+	PetscInt  n = AX[i]->nnp-1; /* last node */
+	A0[i] = 0; if (!w) A1[i] = n;
       }
     }
     {
       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;
       }
     }
   }
       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;
       }
     }
   }
       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;
       }
     }
   }
 /*@
    IGAFormSystem - Form the matrix and vector which represents the
    discretized a(w,u) = L(w).
-   
+
    Collective on IGA/Mat/Vec
 
    Input Parameters:
 +  matA - the matrix obtained from discretization of a(w,u)
 -  vecB - the vector obtained from discretization of L(w)
 
-   Notes: 
+   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.
 {
   IGAUserSystem  System;
   void           *SysCtx;
+  PetscInt       dir,side;
   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);
 #undef  __FUNCT__
 #define __FUNCT__ "IGAFormSystem"
 PetscErrorCode IGAFormSystem(IGA iga,Mat matA,Vec vecB,
-                             IGAUserSystem System,void *ctx)
+			     IGAUserSystem System,void *ctx)
 {
   IGAElement     element;
   IGAPoint       point;
   PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
   IGACheckSetUp(iga,1);
 
-  ierr = MatZeroEntries(matA);CHKERRQ(ierr);
-  ierr = VecZeroEntries(vecB);CHKERRQ(ierr);
-
   /* Element loop */
   ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
   ierr = IGABeginElement(iga,&element);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 (IGABoundaryElementNextPoint(element,point,dir,side)) {
+      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
    communicators as the IGA.
-   
+
    Logically collective on IGA
 
    Input Parameter: