Lisandro Dalcin avatar Lisandro Dalcin committed f0d255a

Rework collocation (periodicity and parallel are broken)

Comments (0)

Files changed (21)

demo/AdvectionDiffusion.c

   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen,dim;
-  IGAPointGetSizes(p,&nen,0,&dim);
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
 
   const PetscReal *N;
   const PetscReal *N1;
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   IGABoundary bnd;
   PetscInt dir,side;
-  for (dir=0; dir<2; dir++) {
+  for (dir=0; dir<3; dir++) {
     for (side=0; side<2; side++) {
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
       ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
   type(IGAPoint), intent(in) :: p
   type(IGAUser),  intent(in) :: ctx
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%nen)
-  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
-  scalar (kind=IGA_SCALAR_KIND ) :: u, grad_u(p%dim)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%neq)
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:), hess_N(:,:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: u, grad_u(p%dim), del2_u
   integer(kind=IGA_INTEGER_KIND) :: a
 
   N      => IGA_Shape0(p)
   grad_N => IGA_Shape1(p)
+  hess_N => IGA_Shape2(p)
 
   u      = IGA_Eval(N,UU)
   grad_u = IGA_Eval(grad_N,UU)
+  del2_u = IGA_Del2(hess_N,UU)
 
   u      = IGA_Value(p,UU) ! just for testing
   grad_u = IGA_Grad(p,UU)  ! just for testing
+  del2_u = IGA_Del2(p,UU)  ! just for testing
 
-  do a = 1, p%nen
-     FF(a) = dot_product(grad_N(:,a),real(grad_u)) - &
+  if (p%neq == 1) then
+     ! collocation
+     FF(1) = - del2_u - ctx%lambda * exp(real(u))
+  else
+     ! galerkin
+     do a = 1, p%nen
+        FF(a) = dot_product(grad_N(:,a),real(grad_u)) - &
              N(a) * ctx%lambda * exp(real(u))
-  end do
+     end do
+  end if
 
   ierr = 0
 end function Function
   type(IGAPoint), intent(in) :: p
   type(IGAUser),  intent(in) :: ctx
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%nen)
-  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
-  scalar (kind=IGA_SCALAR_KIND ) :: u
-  integer(kind=IGA_INTEGER_KIND) :: a, b
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%neq)
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:), hess_N(:,:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: u, grad_u(p%dim), del2_u, del2_N
+  integer(kind=IGA_INTEGER_KIND) :: a, b, i
 
   N      => IGA_Shape0(p)
   grad_N => IGA_Shape1(p)
+  hess_N => IGA_Shape2(p)
 
-  u = IGA_Eval(N,UU)
+  u      = IGA_Eval(N,UU)
+  grad_u = IGA_Eval(grad_N,UU)
+  del2_u = IGA_Del2(hess_N,UU)
 
-  u = IGA_Value(p,UU)
+  u      = IGA_Value(p,UU) ! just for testing
+  grad_u = IGA_Grad(p,UU)  ! just for testing
+  del2_u = IGA_Del2(p,UU)  ! just for testing
 
-  do a = 1, p%nen
-     do b = 1, p%nen
-        JJ(b,a) = dot_product(grad_N(:,a),grad_N(:,b)) - &
-                  N(a) * N(b) * ctx%lambda  * exp(real(u))
+  if (p%neq == 1) then
+     ! collocation
+     do a = 1, p%nen
+        del2_N = 0
+        do i = 1, p%dim
+           del2_N = del2_N + hess_N(i,i,a)
+        end do
+        JJ(a,1) = - del2_N - N(a) * ctx%lambda * exp(real(u))
      end do
-  end do
+  else
+     ! galerkin
+     do a = 1, p%nen
+        do b = 1, p%nen
+           JJ(b,a) = dot_product(grad_N(:,a),grad_N(:,b)) - &
+                N(a) * N(b) * ctx%lambda * exp(real(u))
+        end do
+     end do
+  end if
 
   ierr = 0
 end function Jacobian
   real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%neq)
 
   real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
   scalar (kind=IGA_SCALAR_KIND ) :: v, u, grad_u(p%dim)
   real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%neq)
 
   real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
   scalar (kind=IGA_SCALAR_KIND ) :: v, u, grad_u(p%dim)

demo/CahnHilliard2D.c

   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c_t,c;
   IGAPointFormValue(p,V,&c_t);
   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c_t,c;
   IGAPointFormValue(p,V,&c_t);

demo/CahnHilliard3D.c

   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c,c_t;
   IGAPointFormValue(p,V,&c_t);
   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c_t,c;
   IGAPointFormValue(p,V,&c_t);

demo/L2Projection.c

 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
 PetscErrorCode SystemLaplace(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen,dim;
-  IGAPointGetSizes(p,&nen,0,&dim);
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
 
   const PetscReal *N1;
   IGAPointGetShapeFuns(p,1,&N1);
 PetscErrorCode SystemPoisson(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen,dim;
-  IGAPointGetSizes(p,&nen,0,&dim);
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
 
   const PetscReal *N;
   IGAPointGetShapeFuns(p,0,&N);
 
 #undef  __FUNCT__
 #define __FUNCT__ "SystemCollocation"
-PetscErrorCode SystemCollocation(IGAColPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen,dim;
-  IGAColPointGetSizes(p,&nen,0,&dim);
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
 
   PetscInt Nb[3] = {0,0,0};
-  Nb[0] = p->parent->basis[0]->nnp;
-  Nb[1] = p->parent->basis[1]->nnp;
-  Nb[2] = p->parent->basis[2]->nnp;
+  Nb[0] = p->parent->BD[0]->nnp;
+  Nb[1] = p->parent->BD[1]->nnp;
+  Nb[2] = p->parent->BD[2]->nnp;
 
   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);
+  IGAPointGetBasisFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetBasisFuns(p,1,(const PetscReal**)&N1);
+  IGAPointGetBasisFuns(p,2,(const PetscReal**)&N2);
   
   PetscInt a,i;
   PetscReal n[3] = {0,0,0};
   PetscBool Dirichlet=PETSC_FALSE,Neumann=PETSC_FALSE;
-  for (i=0; i<dim; i++) if (p->ID[i] == 0) Dirichlet=PETSC_TRUE;
+  for (i=0; i<dim; i++) if (p->parent->ID[i] == 0) Dirichlet=PETSC_TRUE;
   if (!Dirichlet) { 
     for (i=0; i<dim; i++) {
-      if (p->ID[i] == Nb[i]-1) {
+      if (p->parent->ID[i] == Nb[i]-1) {
 	Neumann=PETSC_TRUE;
 	n[i] = 1;
 	i=dim;
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   if (Collocation){
-    ierr = IGAColSetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr); 
-    ierr = IGAColComputeSystem(iga,A,b);CHKERRQ(ierr);
+    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr); 
+    ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   }else{
     if (iga->geometry){ 
       ierr = IGASetUserSystem(iga,SystemPoisson,PETSC_NULL);CHKERRQ(ierr); 

demo/PatternFormation.c

   PetscReal tau2  = user->tau2;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
   PetscScalar (*R)[2] = (PetscScalar (*)[2])Re;
 
   PetscScalar uv_t[2],uv_0[2],uv_1[2][2];
   PetscReal tau2  = user->tau2;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
   PetscScalar (*K)[2][nen][2] = (PetscScalar (*)[2][nen][2])Ke;
 
   PetscReal f_u=0,f_v=0;

demo/PhaseFieldCrystal2D.c

   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c_t,c;
   IGAPointFormValue(p,V,&c_t);
   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen;
-  IGAPointGetSizes(p,&nen,0,0);
+  IGAPointGetSizes(p,0,&nen,0);
 
   PetscScalar c_t,c;
   IGAPointFormValue(p,V,&c_t);
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt  nen,dim;
-  IGAPointGetSizes(p,&nen,0,&dim);
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
 runex5b_4:
 	-@${MPIEXEC} -n 4 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
 runex6a_1:
-	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 1 -lambda 3.0
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 1 -lambda 1.0
 runex6a_2:
 	-@${MPIEXEC} -n 2 ./Bratu ${OPTS} -iga_dim 1 -steady false
 runex6b_1:
 	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
 runex6c_4:
 	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
+runex6d_1:
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady -iga_collocation
 
 InputOutput := InputOutput.PETSc runex0a_1 runex0a_2 runex0a_3 runex0a_4 runex0a_8 runex0a.rm InputOutput.rm
 L2Proj1D := L2Proj1D.PETSc runex1a_1 runex1a_4 L2Proj1D.rm
 Poisson1D := Poisson1D.PETSc runex2a_1 runex2a_4 Poisson1D.rm
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
-Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 Bratu.rm
+Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 runex6d_1 Bratu.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 
 /* ---------------------------------------------------------------- */
 
 typedef struct _n_IGAAxis     *IGAAxis;
+typedef struct _n_IGABoundary *IGABoundary;
 typedef struct _n_IGARule     *IGARule;
 typedef struct _n_IGABasis    *IGABasis;
-typedef struct _n_IGABoundary *IGABoundary;
-
 typedef struct _n_IGAElement  *IGAElement;
 typedef struct _n_IGAPoint    *IGAPoint;
 
-typedef struct _n_IGAColPoint *IGAColPoint;
-typedef struct _n_IGAColBasis *IGAColBasis;
-
 /* ---------------------------------------------------------------- */
 
 struct _n_IGAAxis {
 PETSC_EXTERN PetscErrorCode IGABasisDestroy(IGABasis *basis);
 PETSC_EXTERN PetscErrorCode IGABasisReset(IGABasis basis);
 PETSC_EXTERN PetscErrorCode IGABasisReference(IGABasis basis);
-PETSC_EXTERN PetscErrorCode IGABasisInit(IGABasis basis,IGAAxis axis,IGARule rule,PetscInt d);
+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;
 
 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,
   IGAUserSystem     System;
   void              *SysCtx;  
   /**/
-  IGAColUserSystem  ColSystem;
-  void              *ColSysCtx;
-  /**/
   IGAUserFunction   Function;
   void              *FunCtx;
   IGAUserJacobian   Jacobian;
   PetscInt  dof;   /* number of degrees of freedom per node */
   PetscInt  order; /* maximum derivative order */
 
-  IGAAxis  axis[3];
-  IGARule  rule[3];
-  IGABasis basis[3];
+  IGAAxis     axis[3];
+  IGARule     rule[3];
   IGABoundary boundary[3][2];
-  IGAElement  iterator;
 
-  IGAColPoint point_iterator; /* stuff added for collocation */
-  IGAColBasis colbasis[3];
+  IGABasis    elem_basis[3];
+  IGAElement  elem_iterator;
+
+  /* stuff added for collocation */
+  PetscBool   collocation;
+  IGABasis    node_basis[3];
+  IGAElement  node_iterator;
 
   PetscInt  proc_sizes[3];
   PetscInt  proc_ranks[3];
 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 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);
 PETSC_EXTERN PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  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 IGAColSetUserSystem (IGA iga,IGAColUserSystem System,void *SysCtx);
-
 /* ---------------------------------------------------------------- */
 
 struct _n_IGAElement {
   PetscInt index;
   /**/
   PetscInt nqp;
+  PetscInt neq;
   PetscInt nen;
   PetscInt dof;
   PetscInt dim;
   PetscInt nsd;
+  IGABasis *BD;
 
   PetscInt  *mapping;  /*   [nen]      */
 
   IGA      parent;
   IGAPoint iterator;
 
-  PetscInt     nrows,*irows;
-  PetscInt     ncols,*icols;
+  PetscBool   collocation;
+
+  PetscInt    *rowmap;
+  PetscInt    *colmap;
 
   PetscInt     nfix;
   PetscInt    *ifix;
   PetscScalar *wmat[4];
 
 };
+
 PETSC_EXTERN PetscErrorCode IGAElementCreate(IGAElement *element);
 PETSC_EXTERN PetscErrorCode IGAElementDestroy(IGAElement *element);
 PETSC_EXTERN PetscErrorCode IGAElementReset(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementInit(IGAElement element,IGA iga);
 
+PETSC_EXTERN PetscErrorCode IGAGetElement(IGA iga,IGAElement *element);
 PETSC_EXTERN PetscErrorCode IGABeginElement(IGA iga,IGAElement *element);
 PETSC_EXTERN PetscBool      IGANextElement(IGA iga,IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAEndElement(IGA iga,IGAElement *element);
 
-PETSC_EXTERN PetscErrorCode IGAElementInitPoint(IGAElement element,IGAPoint point);
+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 IGAElementEndPoint(IGAElement element,IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAElementBuildShapeFuns(IGAElement element);
 
 PETSC_EXTERN PetscErrorCode IGAElementGetIndex(IGAElement element,PetscInt *index);
-PETSC_EXTERN PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *nen,PetscInt *dof,PetscInt *nqp);
+PETSC_EXTERN PetscErrorCode IGAElementGetCount(IGAElement element,PetscInt *count);
+
+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 *jacobian[],const PetscReal **shapefuns[]);
-
-PETSC_EXTERN PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
+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 IGAElementBuildFix(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementFixValues(IGAElement element,PetscScalar U[]);
+PETSC_EXTERN PetscErrorCode IGAElementFixSystem(IGAElement element,PetscScalar K[],PetscScalar F[]);
 PETSC_EXTERN PetscErrorCode IGAElementFixFunction(IGAElement element,PetscScalar F[]);
 PETSC_EXTERN PetscErrorCode IGAElementFixJacobian(IGAElement element,PetscScalar J[]);
-PETSC_EXTERN PetscErrorCode IGAElementFixSystem(IGAElement element,PetscScalar K[],PetscScalar F[]);
 
 PETSC_EXTERN PetscErrorCode IGAElementAssembleVec(IGAElement element,const PetscScalar F[],Vec vec);
 PETSC_EXTERN PetscErrorCode IGAElementAssembleMat(IGAElement element,const PetscScalar K[],Mat mat);
   PetscInt count;
   PetscInt index;
   /**/
+  PetscInt neq;
   PetscInt nen;
   PetscInt dof;
   PetscInt dim;
                        /*2: [nen][nsd][nsd] */
                        /*3: [nen][nsd][nsd][nsd] */
 
+  IGAElement parent;
+
   PetscInt    nvec;
   PetscScalar *wvec[8];
   PetscInt    nmat;
 PETSC_EXTERN PetscErrorCode IGAPointCreate(IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAPointDestroy(IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAPointReset(IGAPoint point);
+PETSC_EXTERN PetscErrorCode IGAPointInit(IGAPoint point,IGAElement element);
 
 PETSC_EXTERN PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index);
 PETSC_EXTERN PetscErrorCode IGAPointGetCount(IGAPoint point,PetscInt *count);
-PETSC_EXTERN PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim);
+PETSC_EXTERN PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *neq,PetscInt *nen,PetscInt *dof);
+PETSC_EXTERN PetscErrorCode IGAPointGetDims(IGAPoint point,PetscInt *dim,PetscInt *nsd);
 PETSC_EXTERN PetscErrorCode IGAPointGetQuadrature(IGAPoint point,PetscReal *weigth,PetscReal *detJac);
 PETSC_EXTERN PetscErrorCode IGAPointGetBasisFuns(IGAPoint point,PetscInt der,const PetscReal *basisfuns[]);
 PETSC_EXTERN PetscErrorCode IGAPointGetShapeFuns(IGAPoint point,PetscInt der,const PetscReal *shapefuns[]);
 
 /* ---------------------------------------------------------------- */
 
-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 detX;       /*   [1]                  */  
-  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);
-PETSC_EXTERN PetscErrorCode IGAColComputeSystem(IGA iga,Mat A,Vec B);
-
-struct _n_IGAColBasis {
-  PetscInt refct;
-  /**/
-  PetscInt  ncp;      /* number of collocation points */
-  PetscInt  nen;      /* number of local basis functions */
-  PetscInt  p,d;      /* polynomial order, last derivative index */
-
-  PetscInt  *offset;  /* [ncp] basis offset   */
-  PetscReal *detJ;    /* [ncp]                */
-  PetscReal *point;   /* [ncp]                */
-  PetscReal *value;   /* [ncp][nen][d+1]      */
-};
-PETSC_EXTERN PetscErrorCode IGAColBasisCreate(IGAColBasis *basis);
-PETSC_EXTERN PetscErrorCode IGAColBasisDestroy(IGAColBasis *basis);
-PETSC_EXTERN PetscErrorCode IGAColBasisReset(IGAColBasis basis);
-PETSC_EXTERN PetscErrorCode IGAColBasisReference(IGAColBasis basis);
-PETSC_EXTERN PetscErrorCode IGAColBasisInit(IGAColBasis basis,IGAAxis axis,PetscInt d);
-
-/* ---------------------------------------------------------------- */
-
 #ifndef PetscMalloc1
 #define PetscMalloc1(m1,t1,r1) (PetscMalloc((m1)*sizeof(t1),(r1)))
 #endif
 #FFLAGS = ${FFLAGS_STD} -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 petigacol.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
 SOURCEF1 = petigaftn.F90 petigaval.F90
 SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}
   for (i=0; i<3; i++) {
     ierr = IGAAxisCreate(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleCreate(&iga->rule[i]);CHKERRQ(ierr);
-    ierr = IGABasisCreate(&iga->basis[i]);CHKERRQ(ierr);
-    ierr = IGAColBasisCreate(&iga->colbasis[i]);CHKERRQ(ierr); /* collocation */
+    ierr = IGABasisCreate(&iga->elem_basis[i]);CHKERRQ(ierr);
+    ierr = IGABasisCreate(&iga->node_basis[i]);CHKERRQ(ierr);
     ierr = IGABoundaryCreate(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryCreate(&iga->boundary[i][1]);CHKERRQ(ierr);
     iga->proc_sizes[i] = -1;
   }
-  ierr = IGAElementCreate(&iga->iterator);CHKERRQ(ierr);
-  ierr = IGAColPointCreate(&iga->point_iterator);CHKERRQ(ierr);
+  ierr = IGAElementCreate(&iga->elem_iterator);CHKERRQ(ierr);
+  ierr = IGAElementCreate(&iga->node_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
   for (i=0; i<3; i++) {
     ierr = IGAAxisDestroy(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleDestroy(&iga->rule[i]);CHKERRQ(ierr);
-    ierr = IGABasisDestroy(&iga->basis[i]);CHKERRQ(ierr);
-    ierr = IGAColBasisDestroy(&iga->colbasis[i]);CHKERRQ(ierr); /* collocation */
+    ierr = IGABasisDestroy(&iga->elem_basis[i]);CHKERRQ(ierr);
+    ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);
     ierr = IGABoundaryDestroy(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryDestroy(&iga->boundary[i][1]);CHKERRQ(ierr);
   }
-  ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
-  ierr = IGAColPointDestroy(&iga->point_iterator);CHKERRQ(ierr);
+  ierr = IGAElementDestroy(&iga->elem_iterator);CHKERRQ(ierr);
+  ierr = IGAElementDestroy(&iga->node_iterator);CHKERRQ(ierr);
 
   ierr = IGAReset(iga);CHKERRQ(ierr);
   ierr = PetscHeaderDestroy(&iga);CHKERRQ(ierr);
   while (iga->nwork > 0)
     {ierr = VecDestroy(&iga->vwork[--iga->nwork]);CHKERRQ(ierr);}
 
-  ierr = IGAElementReset(iga->iterator);CHKERRQ(ierr);
-  ierr = IGAColPointReset(iga->point_iterator);CHKERRQ(ierr);
+  ierr = IGAElementReset(iga->elem_iterator);CHKERRQ(ierr);
+  ierr = IGAElementReset(iga->node_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
   IGACheckSetUp(iga,1);
   if (i < 0)         SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D must be nonnegative",i);
   if (i >= iga->dim) SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D, but dimension %D",i,iga->dim);
-  *basis = iga->basis[i];
+  *basis = iga->elem_basis[i];
   PetscFunctionReturn(0);
 }
 
         if (nw == 0) w = PETSC_TRUE;
         ierr = IGAAxisSetPeriodic(iga->axis[i],w);CHKERRQ(ierr);
       }
+    
+    /* Collocation */
+    {
+      PetscBool collocation = iga->collocation;
+      ierr = PetscOptionsBool("-iga_collocation","Use collocation","IGASetCollocation",collocation,&collocation,&flg);CHKERRQ(ierr);
+      if (flg) {iga->collocation = collocation;/*ierr = IGASetCollocation(iga,collocation);CHKERRQ(ierr);*/}
+    }
 
     /* Geometry */
     ierr = PetscOptionsString("-iga_geometry","Specify IGA geometry file","IGARead",filename,filename,sizeof(filename),&flg);CHKERRQ(ierr);
     {ierr = IGAAxisSetUp(iga->axis[i]);CHKERRQ(ierr);}
   for (i=iga->dim; i<3; i++) 
     {ierr = IGAAxisReset(iga->axis[i]);CHKERRQ(ierr);}
-
   { /* processor grid and coordinates */
     MPI_Comm    comm = ((PetscObject)iga)->comm;
     PetscMPIInt size,rank;
       geom_gwidth[i] = 1;
     }
   }
-
   PetscFunctionReturn(0);
 }
 
       PetscInt q = iga->axis[i]->p + 1;
       ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);
     }
-    ierr = IGABasisInit(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
+    ierr = IGABasisInitQuadrature(iga->elem_basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
   }
-  ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
+  iga->elem_iterator->collocation = PETSC_FALSE;
+  ierr = IGAElementInit(iga->elem_iterator,iga);CHKERRQ(ierr);
 
   for (i=0; i<3; i++) {
-    ierr = IGAColBasisInit(iga->colbasis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
+    ierr = IGABasisInitCollocation(iga->node_basis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
   }
-  ierr = IGAColPointInit(iga->point_iterator,iga);CHKERRQ(ierr);
-
+  iga->node_iterator->collocation = PETSC_TRUE;
+  ierr = IGAElementInit(iga->node_iterator,iga);CHKERRQ(ierr);
 
   ierr = IGASetUp_View(iga);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAGetElement"
-PetscErrorCode IGAGetElement(IGA iga,IGAElement *element)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidPointer(element,2);
-  *element = iga->iterator;
-  PetscFunctionReturn(0);
-}
-
-#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"
 /*@
    IGASetUserSystem - Set the user callback to form the matrix and vector

src/petigabasis.c

 EXTERN_C_END
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGABasisInit"
-PetscErrorCode IGABasisInit(IGABasis basis,IGAAxis axis,IGARule rule,PetscInt d)
+#define __FUNCT__ "IGABasisInitQuadrature"
+PetscErrorCode IGABasisInitQuadrature(IGABasis basis,IGAAxis axis,IGARule rule,PetscInt d)
 {
   PetscInt       p,nnp;
   const PetscInt *span;
 
   PetscFunctionReturn(0);
 }
+
+EXTERN_C_BEGIN
+extern PetscInt  IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal *U);
+extern PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[]);
+EXTERN_C_END
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABasisInitCollocation"
+PetscErrorCode IGABasisInitCollocation(IGABasis basis,IGAAxis axis,PetscInt d)
+{
+  PetscInt       p,n,nnp;
+  const PetscReal*U;
+  PetscInt       iel,nel;
+  PetscInt       iqp,nqp;
+  PetscInt       nen,ndr;
+  PetscInt       *offset;
+  PetscReal      *detJ;
+  PetscReal      *weight;
+  PetscReal      *point;
+  PetscReal      *value;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(basis,1);
+  PetscValidPointer(axis,2);
+  if (d < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+                      "Derivative order must be grather than zero, got %D",d);
+
+  p = axis->p;
+  n = axis->m - p -1;
+  U = axis->U;
+
+  nel  = axis->nnp;
+  nnp  = axis->nnp;
+  nqp  = 1;
+  nen  = p+1;
+  ndr  = d+1;
+
+  ierr = PetscMalloc1(nel,PetscInt,&offset);CHKERRQ(ierr);
+  ierr = PetscMalloc1(nel,PetscReal,&detJ);CHKERRQ(ierr);
+  ierr = PetscMalloc1(nqp,PetscReal,&weight);CHKERRQ(ierr);
+  ierr = PetscMalloc1(nel*nqp,PetscReal,&point);CHKERRQ(ierr);
+  ierr = PetscMalloc1(nel*nqp*nen*ndr,PetscReal,&value);CHKERRQ(ierr);
+
+  for (iqp=0; iqp<nqp; iqp++) {
+    weight[iqp] = 1.0;
+  }
+  for (iel=0; iel<nel; iel++) {
+    PetscReal u = IGA_Greville(iel,p,U);
+    PetscInt  k = IGA_FindSpan(n,p,u,U);
+    PetscReal *N = &value[iel*nen*ndr];
+    offset[iel] = k-p;
+    point[iel]  = u;
+    detJ[iel]   = U[k+1]-U[k];
+    IGA_DersBasisFuns(k,u,p,d,U,N);
+  }
+
+  ierr = IGABasisReset(basis);CHKERRQ(ierr);
+
+  basis->nel    = nel;
+  basis->nnp    = nnp;
+  basis->nqp    = nqp;
+  basis->nen    = nen;
+  basis->p      = p;
+  basis->d      = d;
+  basis->offset = offset;
+
+  basis->detJ   = detJ;
+  basis->weight = weight;
+  basis->point  = point;
+  basis->value  = value;
+
+  PetscFunctionReturn(0);
+}
+
+PetscInt IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[])
+{
+  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;
+}
+
+PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[])
+{
+  PetscInt j;
+  PetscReal u = 0.0;
+  for(j=0;j<p;j++) u += U[i+j+1];
+  u /= p;
+  return u;
+}

src/petigacol.c

-#include "petiga.h"
-
-extern PetscLogEvent IGA_ColFormSystem;
-
-#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, const 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,const 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)
-{
-  PetscInt i,dim  = point->dim;
-  PetscInt *start = point->start;
-  PetscInt *width = point->width;
-  PetscInt *ID    = point->ID;
-  PetscInt index,coord;
-  IGAColBasis *BD = point->parent->colbasis;
-
-  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]; 
-    point->point[i] = BD[i]->point[ID[i]];
-  }
-  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;
-    IGAColBasis *BD = iga->colbasis;
-    PetscInt *ID = point->ID;
-    
-    PetscInt ia, inen = BD[0]->nen, ioffset = BD[0]->offset[ID[0]];
-    PetscInt ja, jnen = BD[1]->nen, joffset = BD[1]->offset[ID[1]];
-    PetscInt ka, knen = BD[2]->nen, koffset = BD[2]->offset[ID[2]];
-    PetscInt *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(ID,BD,i) 1,BD[i]->nen,BD[i]->d,BD[i]->value+ID[i]*BD[i]->nen*(BD[i]->d+1)
-
-#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;
-  {
-    IGAColBasis *BD  = point->parent->colbasis;
-    PetscInt *ID  = point->ID;
-    PetscReal **N = point->basis;
-    switch (point->dim) {
-    case 3: IGA_BasisFuns_3D(order,point->rational,
-                             point->geometryW,
-                             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,point->rational,
-                             point->geometryW,
-                             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,point->rational,
-                             point->geometryW,
-                             IGA_BasisFuns_ARGS(ID,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->detX;
-    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);
-}
-
-/* The following routines are a slight alteration of the petigaksp.c routines */
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColComputeSystem"
-PetscErrorCode IGAColComputeSystem(IGA iga,Mat matA,Vec vecB)
-{
-  IGAColUserSystem  System;
-  void           *SysCtx;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
-  PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
-  IGACheckSetUp(iga,1);
-  IGACheckUserOp(iga,1,ColSystem);
-  System = iga->userops->ColSystem;
-  SysCtx = iga->userops->ColSysCtx;
-  ierr = IGAColFormSystem(iga,matA,vecB,System,SysCtx);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 = PetscLogEventBegin(IGA_ColFormSystem,iga,matA,vecB,0);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 = PetscLogEventEnd(IGA_ColFormSystem,iga,matA,vecB,0);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);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColSetUserSystem"
-PetscErrorCode IGAColSetUserSystem(IGA iga,IGAColUserSystem System,void *SysCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (System) iga->userops->ColSystem = System;
-  if (SysCtx) iga->userops->ColSysCtx = SysCtx;
-  PetscFunctionReturn(0);
-}
-
-/* The following routines are a slight alteration of the petigabasis.c routines */
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColBasisCreate"
-PetscErrorCode IGAColBasisCreate(IGAColBasis *basis)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(basis,3);
-  ierr = PetscNew(struct _n_IGAColBasis,basis);CHKERRQ(ierr);
-  (*basis)->refct = 1;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColBasisDestroy"
-PetscErrorCode IGAColBasisDestroy(IGAColBasis *_basis)
-{
-  IGAColBasis       basis;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(_basis,1);
-  basis = *_basis; *_basis = 0;
-  if (!basis) PetscFunctionReturn(0);
-  if (--basis->refct > 0) PetscFunctionReturn(0);
-  ierr = IGAColBasisReset(basis);CHKERRQ(ierr);
-  ierr = PetscFree(basis);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColBasisReset"
-PetscErrorCode IGAColBasisReset(IGAColBasis basis)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  if (!basis) PetscFunctionReturn(0);
-  PetscValidPointer(basis,1);
-  basis->ncp = 0;
-  basis->nen = 0;
-  basis->p   = 0;
-  basis->d   = 0;
-  ierr = PetscFree(basis->offset);CHKERRQ(ierr);
-  ierr = PetscFree(basis->detJ);CHKERRQ(ierr);
-  ierr = PetscFree(basis->point);CHKERRQ(ierr);
-  ierr = PetscFree(basis->value);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColBasisReference"
-PetscErrorCode IGAColBasisReference(IGAColBasis basis)
-{
-  PetscFunctionBegin;
-  PetscValidPointer(basis,1);
-  basis->refct++;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAColBasisInit"
-PetscErrorCode IGAColBasisInit(IGAColBasis basis,IGAAxis axis,PetscInt d)
-{
-  PetscInt       m,p;
-  const PetscReal*U;
-  PetscInt       icp,ncp;
-  PetscInt       nen,ndr;
-  PetscInt       *offset;
-  PetscReal      *detJ;
-  PetscReal      *point;
-  PetscReal      *value;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(basis,1);
-  PetscValidPointer(axis,2);
-  if (d < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-                      "Derivative order must be grather than zero, got %D",d);
-
-  m = axis->m;
-  p = axis->p;
-  U = axis->U;
-
-  ncp  = axis->nnp;
-  nen  = p+1;
-  ndr  = d+1;
-
-  ierr = PetscMalloc1(ncp,PetscInt,&offset);CHKERRQ(ierr);
-  ierr = PetscMalloc1(ncp,PetscReal,&detJ);CHKERRQ(ierr);
-  ierr = PetscMalloc1(ncp,PetscReal,&point);CHKERRQ(ierr);
-  ierr = PetscMalloc1(ncp*nen*ndr,PetscReal,&value);CHKERRQ(ierr);
-  
-  for (icp=0; icp<ncp; icp++) {
-    point[icp]   = _Greville(icp,p,m,U); /* XXX Later fix to be a choice */
-    PetscInt k   = _FindSpan(ncp-1,p,point[icp],U);
-    IGA_DersBasisFuns(k,point[icp],p,d,U,&value[icp*nen*ndr]);
-    offset[icp] = k-p;
-    detJ[icp]   = 0.5*(U[k+1]-U[k]);
-  }
-
-  ierr = IGAColBasisReset(basis);CHKERRQ(ierr);
-
-  basis->ncp    = ncp;
-  basis->nen    = nen;
-  basis->p      = p;
-  basis->d      = d;
-  basis->offset = offset;
-
-  basis->detJ   = detJ;
-  basis->point  = point;
-  basis->value  = value;
-
-  PetscFunctionReturn(0);
-}
   ierr = PetscFree(element->shape[2]);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[3]);CHKERRQ(ierr);
 
-  element->nrows = 0; element->irows = PETSC_NULL;
-  element->ncols = 0; element->icols = PETSC_NULL;
   ierr = IGAElementFreeFix(element);CHKERRQ(ierr);
   ierr = IGAElementFreeWork(element);CHKERRQ(ierr);
   ierr = IGAPointReset(element->iterator);CHKERRQ(ierr);
 #define __FUNCT__ "IGAElementInit"
 PetscErrorCode IGAElementInit(IGAElement element,IGA iga)
 {
+  PetscInt *start;
+  PetscInt *width;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   element->dof = iga->dof;
   element->dim = iga->dim;
   element->nsd = iga->nsd ? iga->nsd : iga->dim;
+  if (!element->collocation) {
+    start = iga->elem_start;
+    width = iga->elem_width;
+    element->BD = iga->elem_basis;
+  } else {
+    start = iga->node_lstart;
+    width = iga->node_lwidth;
+    element->BD = iga->node_basis;
+  }
   { /* */
-    IGABasis *BD = iga->basis;
     PetscInt i,dim = element->dim;
-    PetscInt nel=1,nen=1,nqp=1;
+    PetscInt nel=1,nqp=1,nen=1;
     for (i=0; i<dim; i++) {
-      element->start[i] = iga->elem_start[i];
-      element->width[i] = iga->elem_width[i];
+      element->start[i] = start[i];
+      element->width[i] = width[i];
       nel *= element->width[i];
-      nen *= BD[i]->nen;
-      nqp *= BD[i]->nqp;
+      nqp *= element->BD[i]->nqp;
+      nen *= element->BD[i]->nen;
     }
     for (i=dim; i<3; i++) {
       element->start[i] = 0;
     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);
   }
-  { /* */
-    element->nrows = element->nen; element->irows = element->mapping;
-    element->ncols = element->nen; element->icols = element->mapping;
+  if (!element->collocation) {
+    element->neq    = element->nen;
+    element->rowmap = element->mapping;
+    element->colmap = element->mapping;
+  } else {
+    element->neq    = 1;
+    element->rowmap = &element->index;
+    element->colmap = element->mapping;
   }
   { /* */
     PetscInt nen = element->nen;
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->vfix);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->xfix);CHKERRQ(ierr);
   }
-  ierr = IGAElementInitPoint(element,element->iterator);CHKERRQ(ierr);
+  ierr = IGAPointInit(element->iterator,element);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAGetElement"
+PetscErrorCode IGAGetElement(IGA iga,IGAElement *element)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(element,2);
+  if (!iga->collocation)
+    *element = iga->elem_iterator;
+  else
+    *element = iga->node_iterator;
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "IGABeginElement"
 PetscErrorCode IGABeginElement(IGA iga,IGAElement *element)
 {
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidPointer(element,2);
   IGACheckSetUp(iga,1);
-  *element = iga->iterator;
+  ierr = IGAGetElement(iga,element);CHKERRQ(ierr);
   (*element)->index = -1;
   PetscFunctionReturn(0);
 }
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementInitPoint"
-PetscErrorCode IGAElementInitPoint(IGAElement element,IGAPoint point)
+#define __FUNCT__ "IGAElementGetPoint"
+PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point)
 {
-  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidPointer(point,2);
-  ierr = IGAPointReset(point);CHKERRQ(ierr);
-
-  point->nen = element->nen;
-  point->dof = element->dof;
-  point->dim = element->dim;
-  point->nsd = element->nsd;
-
-  point->count = element->nqp;
-  point->index = -1;
-
+  *point = element->iterator;
   PetscFunctionReturn(0);
 }
 
 
   if (element->geometry)
     point->geometry = element->geometryX;
-  else 
+  else
     point->geometry = PETSC_NULL;
   if (element->geometry && dim == nsd) { /* XXX */
     point->detX     = element->detX;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementGetSizes"
-PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *nen,PetscInt *dof,PetscInt *nqp)
+#define __FUNCT__ "IGAElementGetCount"
+PetscErrorCode IGAElementGetCount(IGAElement element,PetscInt *count)
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
-  if (nen) PetscValidIntPointer(nen,2);
-  if (dof) PetscValidIntPointer(dof,3);
-  if (nqp) PetscValidIntPointer(nqp,4);
+  PetscValidIntPointer(count,2);
+  *count = element->count;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementGetSizes"
+PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *neq,PetscInt *nen,PetscInt *dof)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (neq) PetscValidIntPointer(neq,2);
+  if (nen) PetscValidIntPointer(nen,3);
+  if (dof) PetscValidIntPointer(dof,4);
+  if (neq) *neq = element->neq;
   if (nen) *nen = element->nen;
   if (dof) *dof = element->dof;
-  if (nqp) *nqp = element->nqp;
   PetscFunctionReturn(0);
 }
 
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetShapeFuns"
-PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,PetscInt *nen,PetscInt *dim,
-                                      const PetscReal *gradX[],const PetscReal **shapefuns[])
+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 (gradX)     PetscValidPointer(gradX,5);
-  if (shapefuns) PetscValidPointer(shapefuns,6);
+  if (shapefuns) PetscValidPointer(shapefuns,5);
   if (nqp)       *nqp       = element->nqp;
   if (nen)       *nqp       = element->nqp;
   if (dim)       *dim       = element->dim;
-  if (gradX)     *gradX     = element->gradX[0];
   if (shapefuns) *shapefuns = (const PetscReal **)element->shape;
   PetscFunctionReturn(0);
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "IGAElementGetPoint"
-PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point)
-{
-  PetscFunctionBegin;
-  PetscValidPointer(element,1);
-  PetscValidPointer(point,2);
-  *point = element->iterator;
-  PetscFunctionReturn(0);
-}
-
 /* -------------------------------------------------------------------------- */
 
 #undef  __FUNCT__
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   { /* */
     IGA      iga = element->parent;
-    IGABasis *BD = iga->basis;
+    IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
     PetscInt ia, inen = BD[0]->nen, ioffset = BD[0]->offset[ID[0]];
     PetscInt ja, jnen = BD[1]->nen, joffset = BD[1]->offset[ID[1]];
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
-    IGABasis *BD = element->parent->basis;
+    IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
     switch (element->dim) {
     case 3: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   order = element->parent->order;
   {
-    IGABasis *BD  = element->parent->basis;
+    IGABasis *BD  = element->BD;
     PetscInt *ID  = element->ID;
     PetscReal **N = element->basis;
     switch (element->dim) {
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
     size_t MAX_WORK_VEC = sizeof(element->wvec)/sizeof(PetscScalar*);
-    PetscInt m = element->nrows * element->dof;
+    PetscInt m = element->neq * element->dof;
     if (PetscUnlikely(element->nvec >= (PetscInt)MAX_WORK_VEC))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work vectors requested");
     if (PetscUnlikely(!element->wvec[element->nvec])) {
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
     size_t MAX_WORK_MAT = sizeof(element->wmat)/sizeof(PetscScalar*);
-    PetscInt m = element->nrows * element->dof;
-    PetscInt n = element->ncols * element->dof;
+    PetscInt m = element->neq * element->dof;
+    PetscInt n = element->nen * element->dof;
     if (PetscUnlikely(element->nmat >= (PetscInt)MAX_WORK_MAT))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work matrices requested");
     if (PetscUnlikely(!element->wmat[element->nmat])) {
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
-    IGA      iga = element->parent;
-    IGABasis *BD = iga->basis;
+    IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
     PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
     PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
     PetscBool onboundary = PETSC_FALSE;
     PetscInt i,dim = element->dim;
     for (i=0; i<dim; i++) {
-      PetscInt e = BD[i]->nel-1; /* last element */
-      PetscInt n = BD[i]->nnp-1; /* last node    */
-      if (iga->axis[i]->periodic) continue; /* XXX */
-      if (ID[i] == 0) { A0[i] = 0; onboundary = PETSC_TRUE; }
-      if (ID[i] == e) { A1[i] = n; onboundary = PETSC_TRUE; }
+      PetscInt  e = BD[i]->nel-1; /* last element */
+      PetscInt  n = BD[i]->nnp-1; /* last node    */
+      PetscBool b = element->collocation;
+      if (element->parent->axis[i]->periodic) continue; /* XXX */
+      if (ID[i] == 0 || b) { A0[i] = 0; onboundary = PETSC_TRUE; }
+      if (ID[i] == e || b) { A1[i] = n; onboundary = PETSC_TRUE; }
     }
     element->nfix = 0;
     if (onboundary) {
   PetscValidPointer(element,1);
   PetscValidScalarPointer(K,2);
   PetscValidScalarPointer(F,3);
-  {
-    PetscInt M = element->nrows*element->dof;
-    PetscInt N = element->ncols*element->dof;
+  if (!element->collocation) {
+    PetscInt M = element->neq * element->dof;
+    PetscInt N = element->nen * element->dof;
     PetscInt f,nfix=element->nfix;
     for (f=0; f<nfix; f++) {
       PetscInt    k = element->ifix[f];
       PetscInt i,j;
       for (i=0; i<M; i++)
         F[i] -= K[i*N+k] * v;
-      for (i=0; i<M; i++) 
-        K[i*N+k] = 0.0;
-      for (j=0; j<N; j++) 
-        K[k*N+j] = 0.0;
+      for (i=0; i<M; i++) K[i*N+k] = 0.0;
+      for (j=0; j<N; j++) K[k*N+j] = 0.0;
       K[k*N+k] = 1.0;
-      F[k] = v;
+      F[k]     = v;
+    }
+  } else {
+    PetscInt M = element->neq * element->dof;
+    PetscInt N = element->nen * element->dof;
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; 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];
+        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);
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
-  {
+  if (!element->collocation) {
     PetscInt f,nfix=element->nfix;
     for (f=0; f<nfix; f++) {
-      PetscInt    k  = element->ifix[f];
-      PetscScalar u  = element->xfix[f];
-      PetscScalar u0 = element->vfix[f];
-      F[k] = u - u0;
+      PetscInt    k = element->ifix[f];
+      PetscScalar u = element->xfix[f];
+      PetscScalar v = element->vfix[f];
+      F[k] = u - v;
+    }
+  } else {
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      PetscInt k = element->ifix[f];
+      PetscInt a = k / element->dof;
+      if (element->index == element->colmap[a]) {
+        PetscInt    c = k % element->dof;
+        PetscScalar u = element->xfix[f];
+        PetscScalar v = element->vfix[f];
+        F[c] = u - v;
+      }
     }
   }
   PetscFunctionReturn(0);
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(J,2);
-  {
-    PetscInt M = element->nrows*element->dof;
-    PetscInt N = element->ncols*element->dof;
+  if (!element->collocation) {
+    PetscInt M = element->neq * element->dof;
+    PetscInt N = element->nen * element->dof;
     PetscInt f,nfix=element->nfix;
     for (f=0; f<nfix; f++) {
       PetscInt k = element->ifix[f];
       PetscInt i,j;
-      for (i=0; i<M; i++) 
-        J[i*N+k] = 0.0;
-      for (j=0; j<N; j++) 
-        J[k*N+j] = 0.0;
+      for (i=0; i<M; i++) J[i*N+k] = 0.0;
+      for (j=0; j<N; j++) J[k*N+j] = 0.0;
       J[k*N+k] = 1.0;
     }
+  } else {
+    PetscInt M = element->neq * element->dof;
+    PetscInt N = element->nen * element->dof;
+    PetscInt f,nfix=element->nfix;
+    for (f=0; f<nfix; f++) {
+      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;
+      }
+    }
   }
   PetscFunctionReturn(0);
 }
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
   PetscValidHeaderSpecific(vec,VEC_CLASSID,3);
-  mm = element->nrows; ii = element->irows;
+  mm = element->neq; ii = element->rowmap;
   if (element->dof == 1) {
     ierr = VecSetValuesLocal(vec,mm,ii,F,ADD_VALUES);CHKERRQ(ierr);
   } else {
   PetscValidPointer(element,1);
   PetscValidScalarPointer(K,2);
   PetscValidHeaderSpecific(mat,MAT_CLASSID,3);
-  mm = element->nrows; ii = element->irows;
-  nn = element->ncols; jj = element->icols;
+  mm = element->neq; ii = element->rowmap;
+  nn = element->nen; jj = element->colmap;
   if (element->dof == 1) {
     ierr = MatSetValuesLocal(mat,mm,ii,nn,jj,K,ADD_VALUES);CHKERRQ(ierr);
   } else {

src/petigaftn.F90

 
      integer(kind=IGA_INTEGER_KIND) :: count
      integer(kind=IGA_INTEGER_KIND) :: index
+     integer(kind=IGA_INTEGER_KIND) :: neq
      integer(kind=IGA_INTEGER_KIND) :: nen
      integer(kind=IGA_INTEGER_KIND) :: dof
      integer(kind=IGA_INTEGER_KIND) :: dim
   k = i;
   while (U[k]==U[k+1]) k++; /* XXX Using "==" with floating point values ! */
   *first = k - p;
-  /* cmopute index of the rightmost overlapping basis */
+  /* compute index of the rightmost overlapping basis */
   k = i + p + 1;
   while (U[k]==U[k-1]) k--; /* XXX Using "==" with floating point values ! */
   *last = k - 1;
+
+  if (iga->collocation) {
+    PetscInt offset = iga->node_basis[dir]->offset[i];
+    *first = offset;
+    *last  = offset + p;
+  }
 }
 
 PETSC_STATIC_INLINE
   lwidth = iga->node_lwidth;
   for (i=0; i<dim; i++) {
     PetscInt first = lstart[i];
-    PetscInt last  = first + lwidth[i] - 1;
+    PetscInt last  = lstart[i] + lwidth[i] - 1;
     PetscInt gfirst,glast;
     BasisStencil(iga,i,first,&gstart[i],&glast);
     BasisStencil(iga,i,last,&gfirst,&glast);
     ISLocalToGlobalMapping map;
 
     ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-    ierr = IGAElementGetSizes(element,&nen,&dof,0);CHKERRQ(ierr);
+    ierr = IGAElementGetSizes(element,PETSC_NULL,&nen,&dof);CHKERRQ(ierr);
 
     if (dof == 1) {
       ierr = MatGetLocalToGlobalMapping(A,&map,PETSC_NULL);CHKERRQ(ierr);

src/petigapoint.c

 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAPointInit"
+PetscErrorCode IGAPointInit(IGAPoint point,IGAElement element)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  PetscValidPointer(element,2);
+  ierr = IGAPointReset(point);CHKERRQ(ierr);
+  point->parent = element;
+
+  point->neq = element->neq;
+  point->nen = element->nen;
+  point->dof = element->dof;
+  point->dim = element->dim;
+  point->nsd = element->nsd;
+
+  point->count = element->nqp;
+  point->index = -1;
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAPointGetIndex"
 PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index)
 {
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAPointGetSizes"
-PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim)
+PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *neq,PetscInt *nen,PetscInt *dof)
 {
   PetscFunctionBegin;
   PetscValidPointer(point,1);
-  if (nen) PetscValidIntPointer(nen,2);
-  if (dof) PetscValidIntPointer(dof,3);
-  if (dim) PetscValidIntPointer(dim,4);
+  if (neq) PetscValidIntPointer(neq,2);
+  if (nen) PetscValidIntPointer(nen,3);
+  if (dof) PetscValidIntPointer(dof,4);
+  if (neq) *neq = point->neq;
   if (nen) *nen = point->nen;
   if (dof) *dof = point->dof;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointGetDims"
+PetscErrorCode IGAPointGetDims(IGAPoint point,PetscInt *dim,PetscInt *nsd)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(point,1);
+  if (dim) PetscValidIntPointer(dim,2);
+  if (nsd) PetscValidIntPointer(nsd,3);
   if (dim) *dim = point->dim;
+  if (dim) *nsd = point->nsd;
   PetscFunctionReturn(0);
 }
 
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during point loop");
   {
     size_t MAX_WORK_VEC = sizeof(point->wvec)/sizeof(PetscScalar*);
-    PetscInt m = point->nen * point->dof;
+    PetscInt m = point->neq * 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])) {
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during point loop");
   {
     size_t MAX_WORK_MAT = sizeof(point->wmat)/sizeof(PetscScalar*);
-    PetscInt m = point->nen * point->dof;
+    PetscInt m = point->neq * point->dof;
     PetscInt n = point->nen * point->dof;
     if (PetscUnlikely(point->nmat >= (PetscInt)MAX_WORK_MAT))
       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many work matrices requested");