Commits

Nathan Collier committed 4960ad3

first push to implementing collocation methods, lack BCs and mat structure

  • Participants
  • Parent commits 2b0ae81

Comments (0)

Files changed (5)

File demo/Laplace.c

 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "SystemCollocation"
+PetscErrorCode SystemCollocation(IGAColPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen,dim;
+  IGAColPointGetSizes(p,&nen,0,&dim);
+
+  const PetscReal *N0,(*N1)[dim],(*N2)[dim][dim];
+  IGAColPointGetBasisFuns(p,0,(const PetscReal**)&N0);
+  IGAColPointGetBasisFuns(p,1,(const PetscReal**)&N1);
+  IGAColPointGetBasisFuns(p,2,(const PetscReal**)&N2);
+  
+  PetscInt a,i;
+  for (a=0; a<nen; a++) 
+    for (i=0; i<dim; i++) {
+      K[a] += -N2[a][i][i];
+    }
+      
+  return 0;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "ErrorLaplace"
 PetscErrorCode ErrorLaplace(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
 {
   PetscInt  dim = 3; 
   PetscBool print_error = PETSC_FALSE; 
   PetscBool draw = PETSC_FALSE; 
+  PetscBool Collocation = PETSC_FALSE;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace Options","IGA");CHKERRQ(ierr); 
   ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr); 
   ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr); 
   ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr); 
+  ierr = PetscOptionsBool("-collocation","Enable to use collocation",__FILE__,Collocation,&Collocation,PETSC_NULL);CHKERRQ(ierr); 
   ierr = PetscOptionsEnd();CHKERRQ(ierr); 
 
   // Initialize the discretization
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGAView(iga,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-
   // Set boundary conditions
-
+  
   if (iga->geometry) {
     for(i=0; i<dim; i++) {
       IGABoundary bnd;
       ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
     }
   }
-  
+    
   // Assemble
 
   Mat A;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  if (iga->geometry){
-    ierr = IGAFormSystem(iga,A,b,SystemPoisson,PETSC_NULL);CHKERRQ(ierr);
-  }else{
-    ierr = IGAFormSystem(iga,A,b,SystemLaplace,PETSC_NULL);CHKERRQ(ierr);
-  }
-  ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+  if ( iga->geometry && !Collocation){ ierr = IGAFormSystem(iga,A,b,SystemPoisson,PETSC_NULL);CHKERRQ(ierr); }
+  if (!iga->geometry && !Collocation){ ierr = IGAFormSystem(iga,A,b,SystemLaplace,PETSC_NULL);CHKERRQ(ierr); }
+  if (!iga->geometry &&  Collocation){ ierr = IGAColFormSystem(iga,A,b,SystemCollocation,PETSC_NULL);CHKERRQ(ierr); }
 
   // Solve
 

File include/petiga.h

 
 typedef struct _n_IGAElement  *IGAElement;
 typedef struct _n_IGAPoint    *IGAPoint;
+typedef struct _n_IGAColPoint *IGAColPoint;
 
 /* ---------------------------------------------------------------- */
 
 
 typedef PetscErrorCode (*IGAUserScalar)    (IGAPoint point,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx);
 typedef PetscErrorCode (*IGAUserSystem)    (IGAPoint point,PetscScalar *K,PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAColUserSystem) (IGAColPoint point,PetscScalar *K,PetscScalar *F,void *ctx);
 typedef PetscErrorCode (*IGAUserFunction)  (IGAPoint point,const PetscScalar *U,PetscScalar *F,void *ctx);
 typedef PetscErrorCode (*IGAUserJacobian)  (IGAPoint point,const PetscScalar *U,PetscScalar *J,void *ctx);
 typedef PetscErrorCode (*IGAUserIFunction) (IGAPoint point,PetscReal dt,
 struct _IGAUserOps {
   /**/
   IGAUserSystem     System;
-  void              *SysCtx;
+  void              *SysCtx;  
+  /**/
+  IGAColUserSystem  ColSystem;
+  void              *ColSysCtx;
   /**/
   IGAUserFunction   Function;
   void              *FunCtx;
   IGABasis basis[3];
   IGABoundary boundary[3][2];
   IGAElement  iterator;
+  IGAColPoint point_iterator;
 
   PetscInt  proc_sizes[3];
   PetscInt  proc_ranks[3];
 PETSC_EXTERN PetscErrorCode IGARestoreLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 
 PETSC_EXTERN PetscErrorCode IGAGetElement(IGA iga,IGAElement *element);
+PETSC_EXTERN PetscErrorCode IGAGetColPoint(IGA iga,IGAColPoint *point);
 
 PETSC_EXTERN PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *SysCtx);
 PETSC_EXTERN PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *FunCtx);
   PetscReal *weight;   /*   [1]   */
   PetscReal *detJac;   /*   [1]   */
 
-  PetscReal *point;    /*   [dim] */
-  PetscReal *scale;    /*   [dim] */
-  PetscReal *basis[4]; /*0: [nen] */
-                       /*1: [nen][dim] */
-                       /*2: [nen][dim][dim] */
-                       /*3: [nen][dim][dim][dim] */
-
-  PetscReal *gradX[2]; /*0: [nsd][dim] */
-                       /*1: [dim][nsd] */
-  PetscReal *shape[4]; /*0: [nen]  */
-                       /*1: [nen][nsd] */
-                       /*2: [nen][nsd][nsd] */
-                       /*3: [nen][nsd][nsd][nsd] */
+  PetscReal *point;    /*   [1][dim] */
+  PetscReal *scale;    /*   [1][dim] */
+  
+  PetscReal *basis[4]; /*0: [1][nen] */
+                       /*1: [1][nen][dim] */
+                       /*2: [1][nen][dim][dim] */
+                       /*3: [1][nen][dim][dim][dim] */
+  PetscReal *gradX[2]; /*0: [1][nsd][dim] */
+                       /*1: [1][dim][nsd] */
+  PetscReal *shape[4]; /*0: [1][nen]  */
+                       /*1: [1][nen][nsd] */
+                       /*2: [1][nen][nsd][nsd] */
+                       /*3: [1][nen][nsd][nsd][nsd] */
 
   IGAElement  parent;
 
 
 /* ---------------------------------------------------------------- */
 
+struct _n_IGAColPoint {
+  PetscInt refct;
+  /**/
+  PetscInt count;
+  PetscInt index;  
+  PetscInt start[3];
+  PetscInt width[3];
+  PetscInt ID[3]; 
+  PetscInt span[3];
+  
+  PetscInt nen;
+  PetscInt dof;
+  PetscInt dim;
+  PetscInt nsd;
+  
+  PetscInt *mapping;    /*  [nen] */
+  
+  PetscBool geometry;
+  PetscBool rational;
+  PetscReal *geometryX; /*   [nen][nsd] */
+  PetscReal *geometryW; /*   [nen]      */
+
+  PetscReal detJac;  
+  PetscReal *point;     /*   [dim]                */
+  PetscReal *scale;     /*   [dim]                */
+  
+  PetscReal *basis1d[3];
+
+  PetscReal *basis[4];  /*0: [nen]                */
+                        /*1: [nen][dim]           */
+                        /*2: [nen][dim][dim]      */
+                        /*3: [nen][dim][dim][dim] */
+
+  PetscReal *gradX[2];  /*0: [nsd][dim]           */
+                        /*1: [dim][nsd]           */
+  PetscReal *shape[4];  /*0: [nen]                */
+                        /*1: [nen][nsd]           */
+                        /*2: [nen][nsd][nsd]      */
+                        /*3: [nen][nsd][nsd][nsd] */
+
+  PetscInt    nvec;
+  PetscScalar *wvec[8];
+  PetscInt    nmat;
+  PetscScalar *wmat[4];
+
+  IGA      parent;
+};
+
+PETSC_EXTERN PetscErrorCode IGAColPointCreate(IGAColPoint *point);
+PETSC_EXTERN PetscErrorCode IGAColPointDestroy(IGAColPoint *point);
+PETSC_EXTERN PetscErrorCode IGAColPointReset(IGAColPoint point);
+PETSC_EXTERN PetscErrorCode IGAColPointInit(IGAColPoint point,IGA iga);
+
+PETSC_EXTERN PetscErrorCode IGAColPointBegin(IGAColPoint point);
+PETSC_EXTERN PetscBool      IGAColPointNext(IGAColPoint point);
+PETSC_EXTERN PetscErrorCode IGAColPointEnd(IGAColPoint point);
+
+PETSC_EXTERN PetscErrorCode IGAColPointBuildMapping(IGAColPoint point);
+PETSC_EXTERN PetscErrorCode IGAColPointBuildGeometry(IGAColPoint point);
+PETSC_EXTERN PetscErrorCode IGAColPointBuildShapeFuns(IGAColPoint point);
+
+PETSC_EXTERN PetscErrorCode IGAColPointGetIndex(IGAColPoint point,PetscInt *index);
+PETSC_EXTERN PetscErrorCode IGAColPointGetSizes(IGAColPoint point,PetscInt *nen,PetscInt *dof,PetscInt *nqp);
+PETSC_EXTERN PetscErrorCode IGAColPointGetSizes(IGAColPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim);
+PETSC_EXTERN PetscErrorCode IGAColPointGetMapping(IGAColPoint point,PetscInt *nen,const PetscInt *mapping[]);
+PETSC_EXTERN PetscErrorCode IGAColPointGetShapeFuns(IGAColPoint point,PetscInt der,const PetscReal *shapefuns[]);
+PETSC_EXTERN PetscErrorCode IGAColPointGetBasisFuns(IGAColPoint point,PetscInt der,const PetscReal *basisfuns[]);
+
+PETSC_EXTERN PetscErrorCode IGAColPointGetWorkVec(IGAColPoint point,PetscScalar *V[]);
+PETSC_EXTERN PetscErrorCode IGAColPointGetWorkMat(IGAColPoint point,PetscScalar *M[]);
+PETSC_EXTERN PetscErrorCode IGAColPointGetValues(IGAColPoint point,const PetscScalar U[],PetscScalar u[]);
+
+PETSC_EXTERN PetscErrorCode IGAColPointAssembleVec(IGAColPoint point,const PetscScalar F[],Vec vec);
+PETSC_EXTERN PetscErrorCode IGAColPointAssembleMat(IGAColPoint point,const PetscScalar K[],Mat mat);
+PETSC_EXTERN PetscErrorCode IGAColFormSystem(IGA iga,Mat matA,Vec vecB,IGAColUserSystem System,void *ctx);
+
+/* ---------------------------------------------------------------- */
+
 #ifndef PetscMalloc1
 #define PetscMalloc1(m1,t1,r1) (PetscMalloc((m1)*sizeof(t1),(r1)))
 #endif

File src/makefile

 #FFLAGS = -g3 -pedantic -Wall -Wextra -fcheck=all
 
 SOURCEH  = ../include/petiga.h petigapc.h petigagrid.h petigapart.h
-SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigaio.c petigagrid.c petigapart.c
+SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigaio.c petigagrid.c petigapart.c petigacol.c
 SOURCEF1 = petigaftn.F90
 SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90 petigaint.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}

File src/petiga.c

     iga->proc_sizes[i] = -1;
   }
   ierr = IGAElementCreate(&iga->iterator);CHKERRQ(ierr);
+  ierr = IGAColPointCreate(&iga->point_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
     ierr = IGABoundaryDestroy(&iga->boundary[i][1]);CHKERRQ(ierr);
   }
   ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
+  ierr = IGAColPointDestroy(&iga->point_iterator);CHKERRQ(ierr);
 
   ierr = IGAReset(iga);CHKERRQ(ierr);
   ierr = PetscHeaderDestroy(&iga);CHKERRQ(ierr);
     {ierr = VecDestroy(&iga->vwork[--iga->nwork]);CHKERRQ(ierr);}
 
   ierr = IGAElementReset(iga->iterator);CHKERRQ(ierr);
+  ierr = IGAColPointReset(iga->point_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
     ierr = IGABasisInit(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
   }
   ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
+  ierr = IGAColPointInit(iga->point_iterator,iga);CHKERRQ(ierr);
 
 
   ierr = IGASetUp_View(iga);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAGetColPoint"
+PetscErrorCode IGAGetColPoint(IGA iga,IGAColPoint *point)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(point,2);
+  *point = iga->point_iterator;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGASetUserSystem"
 PetscErrorCode IGASetUserSystem(IGA iga,IGAUserSystem System,void *SysCtx)
 {

File src/petigacol.c

+#include "petiga.h"
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointCreate"
+PetscErrorCode IGAColPointCreate(IGAColPoint *_point)
+{
+  IGAColPoint point;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(_point,1);
+  ierr = PetscNew(struct _n_IGAColPoint,_point);CHKERRQ(ierr);
+  point = *_point;
+  point->refct =  1;
+  point->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointDestroy"
+PetscErrorCode IGAColPointDestroy(IGAColPoint *_point)
+{
+  IGAColPoint     point;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(_point,1);
+  point = *_point; *_point = 0;
+  if (!point) PetscFunctionReturn(0);
+  if (--point->refct > 0) PetscFunctionReturn(0);
+  ierr = IGAColPointReset(point);CHKERRQ(ierr);
+  ierr = PetscFree(point);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointFreeWork"
+static
+PetscErrorCode IGAColPointFreeWork(IGAColPoint point)
+{
+  size_t i;
+  size_t MAX_WORK_VEC;
+  size_t MAX_WORK_MAT;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  MAX_WORK_VEC = sizeof(point->wvec)/sizeof(PetscScalar*);
+  for (i=0; i<MAX_WORK_VEC; i++)
+    {ierr = PetscFree(point->wvec[i]);CHKERRQ(ierr);}
+  point->nvec = 0;
+  MAX_WORK_MAT = sizeof(point->wmat)/sizeof(PetscScalar*);
+  for (i=0; i<MAX_WORK_MAT; i++)
+    {ierr = PetscFree(point->wmat[i]);CHKERRQ(ierr);}
+  point->nmat = 0;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointReset"
+PetscErrorCode IGAColPointReset(IGAColPoint point)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  if (!point) PetscFunctionReturn(0);
+  PetscValidPointer(point,1);
+  point->index = -1;
+  ierr = IGAColPointFreeWork(point);CHKERRQ(ierr);
+
+  ierr = PetscFree(point->mapping);CHKERRQ(ierr);
+  ierr = PetscFree(point->geometryX);CHKERRQ(ierr);
+  ierr = PetscFree(point->geometryW);CHKERRQ(ierr);
+
+  ierr = PetscFree(point->point);CHKERRQ(ierr);
+  ierr = PetscFree(point->scale);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis1d[0]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis1d[1]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis1d[2]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis[0]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis[1]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis[2]);CHKERRQ(ierr);
+  ierr = PetscFree(point->basis[3]);CHKERRQ(ierr);
+
+  ierr = PetscFree(point->gradX[0]);CHKERRQ(ierr);
+  ierr = PetscFree(point->gradX[1]);CHKERRQ(ierr);
+  ierr = PetscFree(point->shape[0]);CHKERRQ(ierr);
+  ierr = PetscFree(point->shape[1]);CHKERRQ(ierr);
+  ierr = PetscFree(point->shape[2]);CHKERRQ(ierr);
+  ierr = PetscFree(point->shape[3]);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointInit"
+PetscErrorCode IGAColPointInit(IGAColPoint point,IGA iga)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
+  if (PetscUnlikely(!iga->setup)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetUp() first");
+  ierr = IGAColPointReset(point);CHKERRQ(ierr);
+  point->parent = iga;
+
+  point->dof = iga->dof;
+  point->dim = iga->dim;
+  point->nsd = iga->nsd ? iga->nsd : iga->dim;
+
+  IGABasis *BD = iga->basis;
+  { /* */
+    PetscInt i,dim = point->dim;
+    PetscInt npts=1,nen=1;
+    for (i=0; i<dim; i++) {
+      point->start[i] = iga->node_lstart[i];
+      point->width[i] = iga->node_lwidth[i];
+      npts *= iga->node_lwidth[i];
+      nen  *= BD[i]->nen;
+    }
+    point->index = -1;
+    point->count = npts;
+    point->nen   = nen;
+  }
+  { /* */
+    PetscInt dim = point->dim;
+    PetscInt nsd = point->nsd;
+    PetscInt nen = point->nen;
+    PetscInt nqp = 1; // looks like the element code, but just for 1 point
+
+    ierr = PetscMalloc1(nen,PetscInt,&point->mapping);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*nsd,PetscReal,&point->geometryX);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen    ,PetscReal,&point->geometryW);CHKERRQ(ierr);    
+
+    ierr = PetscMalloc1(nqp*dim,PetscReal,&point->point);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*dim,PetscReal,&point->scale);CHKERRQ(ierr);
+
+    ierr = PetscMalloc1(BD[0]->nen*(iga->order+1),PetscReal,&point->basis1d[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(BD[1]->nen*(iga->order+1),PetscReal,&point->basis1d[1]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(BD[2]->nen*(iga->order+1),PetscReal,&point->basis1d[2]);CHKERRQ(ierr);
+
+    ierr = PetscMalloc1(nqp*nen,PetscReal,&point->basis[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim,PetscReal,&point->basis[1]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&point->basis[2]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&point->basis[3]);CHKERRQ(ierr);
+
+    ierr = PetscMalloc1(nqp*dim*dim,PetscReal,&point->gradX[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*dim*dim,PetscReal,&point->gradX[1]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen,PetscReal,&point->shape[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim,PetscReal,&point->shape[1]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&point->shape[2]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&point->shape[3]);CHKERRQ(ierr);
+
+    ierr = PetscMemzero(point->basis1d[0],BD[0]->nen*(iga->order+1)*sizeof(PetscReal));CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis1d[1],BD[1]->nen*(iga->order+1)*sizeof(PetscReal));CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis1d[2],BD[2]->nen*(iga->order+1)*sizeof(PetscReal));CHKERRQ(ierr);
+
+    ierr = PetscMemzero(point->point,   sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->scale,   sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis[0],sizeof(PetscReal)*nqp*nen);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis[1],sizeof(PetscReal)*nqp*nen*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->basis[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
+
+    ierr = PetscMemzero(point->gradX[0],sizeof(PetscReal)*nqp*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->gradX[1],sizeof(PetscReal)*nqp*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->shape[0],sizeof(PetscReal)*nqp*nen);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->shape[1],sizeof(PetscReal)*nqp*nen*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->shape[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(point->shape[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointBegin"
+PetscErrorCode IGAColPointBegin(IGAColPoint point)
+{
+  IGA            iga;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  iga = point->parent;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
+  if (PetscUnlikely(!iga->setup)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetUp() first");
+  point->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "_FindSpan"
+PetscInt _FindSpan(PetscInt n,PetscInt p,PetscReal u, PetscReal *U)
+{
+  /* n is the index of the last basis */
+  PetscInt low,high,span;
+  if(u >= U[n+1]) return n;
+  if(u <= U[p]) return p;
+  low  = p;
+  high = n+1;
+  span = (high+low)/2;
+  while(u < U[span] || u >= U[span+1]){
+    if(u < U[span])
+      high = span;
+    else
+      low = span;
+    span = (high+low)/2;
+  }
+  return span;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "_Greville"
+PetscReal _Greville(PetscInt i,PetscInt p,PetscInt m,PetscReal *U)
+{
+  PetscInt j;
+  PetscReal X = 0.0;
+  for(j=0;j<p;j++) X += U[i+j+1];
+  X /= p;
+  return X;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointNext"
+PetscBool IGAColPointNext(IGAColPoint point)
+{
+  IGA      iga = point->parent;
+  PetscInt i,dim  = point->dim;
+  PetscInt *start = point->start;
+  PetscInt *width = point->width;
+  PetscInt *ID    = point->ID;
+  PetscInt *span  = point->span;
+  PetscReal *pnt  = point->point;
+  PetscInt index,coord;
+ 
+  point->nvec = 0;
+  point->nmat = 0;
+
+  index = ++point->index;
+  if (PetscUnlikely(index >= point->count)) {
+    point->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]; 
+    pnt[i]  = _Greville(ID[i],iga->axis[i]->p,iga->axis[i]->m,iga->axis[i]->U);
+    span[i] = _FindSpan(iga->axis[i]->m-iga->axis[i]->p-1,iga->axis[i]->p,pnt[i],iga->axis[i]->U);
+  }
+  for (i=dim; i<3; i++) {
+    span[i] = 0;
+  }
+  //printf("ID: {%d,%d,%d}  pnt: {%.2f,%.2f,%.2f}  span: {%d,%d,%d}\n",ID[0],ID[1],ID[2],pnt[0],pnt[1],pnt[2],span[0],span[1],span[2]);
+  IGAColPointBuildMapping(point);
+  IGAColPointBuildGeometry(point);
+  IGAColPointBuildShapeFuns(point);
+  return PETSC_TRUE;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointEnd"
+PetscErrorCode IGAColPointEnd(IGAColPoint point)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  point->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointBuildMapping"
+PetscErrorCode IGAColPointBuildMapping(IGAColPoint point)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  if (PetscUnlikely(point->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during collocation point loop");
+  { /* */
+    IGA      iga = point->parent;
+    IGABasis *BD = iga->basis;
+    PetscInt *span = point->span;
+    PetscInt ia, inen = BD[0]->nen, ioffset = span[0]-iga->axis[0]->p;
+    PetscInt ja, jnen = BD[1]->nen, joffset = span[1]-iga->axis[1]->p;
+    PetscInt ka, knen = BD[2]->nen, koffset = span[2]-iga->axis[2]->p;
+    PetscInt *start = iga->node_gstart, *width = iga->node_gwidth;
+    PetscInt istart = start[0]/*istride = 1*/;
+    PetscInt jstart = start[1], jstride = width[0];
+    PetscInt kstart = start[2], kstride = width[0]*width[1];
+    PetscInt a=0, *mapping = point->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;
+        }
+      }
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointBuildGeometry"
+PetscErrorCode IGAColPointBuildGeometry(IGAColPoint point)
+{
+  IGA iga;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  if (PetscUnlikely(point->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during collocation point loop");
+  iga = point->parent;
+  point->geometry = iga->geometry;
+  point->rational = iga->rational;
+  if (point->geometry || point->rational) {
+    const PetscInt  *map = point->mapping;
+    const PetscReal *arrayX = iga->geometryX;
+    const PetscReal *arrayW = iga->geometryW;
+    PetscReal *X = point->geometryX;
+    PetscReal *W = point->geometryW;
+    PetscInt a,nen = point->nen;
+    PetscInt i,nsd = point->nsd;
+    if (point->geometry)
+      for (a=0; a<nen; a++)
+        for (i=0; i<nsd; i++)
+          X[i + a*nsd] = arrayX[i + map[a]*nsd];
+    if (point->rational)
+      for (a=0; a<nen; a++)
+        W[a] = arrayW[map[a]];
+  }
+  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
+
+EXTERN_C_BEGIN
+extern void IGA_DersBasisFuns(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
+EXTERN_C_END
+
+#define IGA_BasisFuns_ARGS(BD,i) 1,BD[i]->nen,BD[i]->d,point->basis1d[i]
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointBuildShapeFuns"
+PetscErrorCode IGAColPointBuildShapeFuns(IGAColPoint point)
+{
+  PetscInt order;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  if (PetscUnlikely(point->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during collocation point loop");
+  order = point->parent->order;
+  {
+    IGABasis *BD  = point->parent->basis;
+    PetscReal **N = point->basis;
+    switch (point->dim) {
+    case 3: 
+      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
+      IGA_DersBasisFuns(point->span[1],point->point[1],point->parent->axis[1]->p,point->parent->order,point->parent->axis[1]->U,point->basis1d[1]);
+      IGA_DersBasisFuns(point->span[2],point->point[2],point->parent->axis[2]->p,point->parent->order,point->parent->axis[2]->U,point->basis1d[2]);
+      IGA_BasisFuns_3D(order,point->rational,
+		       point->geometryW,
+		       IGA_BasisFuns_ARGS(BD,0),
+		       IGA_BasisFuns_ARGS(BD,1),
+		       IGA_BasisFuns_ARGS(BD,2),
+		       N[0],N[1],N[2],N[3]); break;
+    case 2: 
+      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
+      IGA_DersBasisFuns(point->span[1],point->point[1],point->parent->axis[1]->p,point->parent->order,point->parent->axis[1]->U,point->basis1d[1]);
+      IGA_BasisFuns_2D(order,point->rational,
+		       point->geometryW,
+		       IGA_BasisFuns_ARGS(BD,0),
+		       IGA_BasisFuns_ARGS(BD,1),
+		       N[0],N[1],N[2],N[3]); break;
+    case 1: 
+      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
+      IGA_BasisFuns_1D(order,point->rational,
+		       point->geometryW,
+		       IGA_BasisFuns_ARGS(BD,0),
+		       N[0],N[1],N[2],N[3]); break;
+    }
+  }
+  if (point->dim == point->nsd) /* XXX */
+  if (point->geometry) {
+    PetscReal **M = point->basis;
+    PetscReal **N = point->shape;
+    PetscReal dX  = point->detJac;
+    PetscReal **gX = point->gradX;
+    switch (point->dim) {
+    case 3: IGA_ShapeFuns_3D(order,
+                             1,point->nen,
+                             point->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,
+                             1,point->nen,
+                             point->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,
+                             1,point->nen,
+                             point->geometryX,
+                             M[0],M[1],M[2],M[3],
+                             N[0],N[1],N[2],N[3],
+                             &dX,gX[0],gX[1]); break;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetIndex"
+PetscErrorCode IGAColPointGetIndex(IGAColPoint point,PetscInt *index)
+{
+  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Function not implemented");
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetSizes"
+PetscErrorCode IGAColPointGetSizes(IGAColPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  if (nen) PetscValidIntPointer(nen,2);
+  if (dof) PetscValidIntPointer(dof,3);
+  if (dim) PetscValidIntPointer(dim,4);
+  if (nen) *nen = point->nen;
+  if (dof) *dof = point->dof;
+  if (dim) *dim = point->dim;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetMapping"
+PetscErrorCode IGAColPointGetMapping(IGAColPoint point,PetscInt *nen,const PetscInt *mapping[])
+{
+  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Function not implemented");
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetShapeFuns"
+PetscErrorCode IGAColPointGetShapeFuns(IGAColPoint point,PetscInt der,const PetscReal *shapefuns[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidPointer(shapefuns,2);
+  if (PetscUnlikely(der < 0 || der >= (PetscInt)(sizeof(point->shape)/sizeof(PetscReal*))))
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+            "Requested derivative must be in range [0,%d], got %D",
+            (int)(sizeof(point->shape)/sizeof(PetscReal*)-1),der);
+  *shapefuns = point->shape[der];
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetBasisFuns"
+PetscErrorCode IGAColPointGetBasisFuns(IGAColPoint point,PetscInt der,const PetscReal *basisfuns[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidPointer(basisfuns,2);
+  if (PetscUnlikely(der < 0 || der >= (PetscInt)(sizeof(point->basis)/sizeof(PetscReal*))))
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+            "Requested derivative must be in range [0,%d], got %D",
+            (int)(sizeof(point->basis)/sizeof(PetscReal*)-1),der);
+  *basisfuns = point->basis[der];
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetWorkVec"
+PetscErrorCode IGAColPointGetWorkVec(IGAColPoint point,PetscScalar *V[])
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidPointer(V,2);
+  if (PetscUnlikely(point->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during colocation point loop");
+  {
+    size_t MAX_WORK_VEC = sizeof(point->wvec)/sizeof(PetscScalar*);
+    PetscInt n = point->dof;
+    if (PetscUnlikely(point->nvec >= (PetscInt)MAX_WORK_VEC))
+      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work vectors requested");
+    if (PetscUnlikely(!point->wvec[point->nvec])) {
+      ierr = PetscMalloc1(n,PetscScalar,&point->wvec[point->nvec]);CHKERRQ(ierr);
+    }
+    *V = point->wvec[point->nvec++];
+    ierr = PetscMemzero(*V,n*sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetWorkMat"
+PetscErrorCode IGAColPointGetWorkMat(IGAColPoint point,PetscScalar *M[])
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidPointer(M,2);
+  if (PetscUnlikely(point->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during colocation point loop");
+  {
+    size_t MAX_WORK_MAT = sizeof(point->wmat)/sizeof(PetscScalar*);
+    PetscInt n = point->nen * point->dof;
+    if (PetscUnlikely(point->nmat >= (PetscInt)MAX_WORK_MAT))
+      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work matrices requested");
+    if (PetscUnlikely(!point->wmat[point->nmat])) {
+      ierr = PetscMalloc1(n,PetscScalar,&point->wmat[point->nmat]);CHKERRQ(ierr);
+    }
+    *M = point->wmat[point->nmat++];
+    ierr = PetscMemzero(*M,n*sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointGetValues"
+PetscErrorCode IGAColPointGetValues(IGAColPoint point,const PetscScalar U[],PetscScalar u[])
+{
+  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Function not implemented");
+  PetscFunctionReturn(0);
+}
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointAssembleVec"
+PetscErrorCode IGAColPointAssembleVec(IGAColPoint point,const PetscScalar F[],Vec vec)
+{
+  PetscInt       index;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidScalarPointer(F,2);
+  PetscValidHeaderSpecific(vec,VEC_CLASSID,3);
+  index = point->index;
+  if (point->dof == 1) {
+    ierr = VecSetValuesLocal(vec,1,&index,F,ADD_VALUES);CHKERRQ(ierr);
+  } else {
+    ierr = VecSetValuesBlockedLocal(vec,1,&index,F,ADD_VALUES);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColPointAssembleMat"
+PetscErrorCode IGAColPointAssembleMat(IGAColPoint point,const PetscScalar K[],Mat mat)
+{
+  PetscInt       nn,index,*ii;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidScalarPointer(K,2);
+  PetscValidHeaderSpecific(mat,MAT_CLASSID,3);
+  nn = point->nen;
+  ii = point->mapping;
+  index = point->index;
+  if (point->dof == 1) {
+    ierr = MatSetValuesLocal(mat,1,&index,nn,ii,K,ADD_VALUES);CHKERRQ(ierr);
+  } else {
+    ierr = MatSetValuesBlockedLocal(mat,1,&index,nn,ii,K,ADD_VALUES);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColFormSystem"
+PetscErrorCode IGAColFormSystem(IGA iga,Mat matA,Vec vecB,IGAColUserSystem System,void *ctx)
+{
+  IGAColPoint     point;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
+  PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
+  IGACheckSetUp(iga,1);
+  
+  ierr = MatZeroEntries(matA);CHKERRQ(ierr);
+  ierr = VecZeroEntries(vecB);CHKERRQ(ierr);
+  ierr = IGAGetColPoint(iga,&point);CHKERRQ(ierr);
+  ierr = IGAColPointBegin(point);CHKERRQ(ierr);
+  while (IGAColPointNext(point)) {
+    PetscScalar *K, *F;
+    ierr = IGAColPointGetWorkMat(point,&K);CHKERRQ(ierr);
+    ierr = IGAColPointGetWorkVec(point,&F);CHKERRQ(ierr);
+    ierr = System(point,K,F,ctx);CHKERRQ(ierr);
+    ierr = IGAColPointAssembleMat(point,K,matA);CHKERRQ(ierr);
+    ierr = IGAColPointAssembleVec(point,F,vecB);CHKERRQ(ierr);
+  }
+  ierr = IGAColPointEnd(point);CHKERRQ(ierr);
+
+  ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd  (matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(vecB);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd  (vecB);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}