Commits

Lisandro Dalcin committed 8e686a0

Implement geometrical mapping

Comments (0)

Files changed (11)

   IGARule  rule[3];
   IGABasis basis[3];
   IGABoundary boundary[3][2];
+  IGAElement iterator;
+
+  const PetscScalar *geometry;
+  PetscBool          rational;
+  Vec                vec_geom;
+  DM                 dm_geom;
 
   PetscInt proc_rank[3];
   PetscInt proc_sizes[3];
   PetscInt elem_sizes[3];
   PetscInt elem_start[3];
   PetscInt elem_width[3];
-
-  DM dm_geom;
-  PetscBool rational;
-  PetscBool geometry;
-  DM dm_dof;
-
-  IGAElement iterator;
+  DM       dm_dof;
 };
 
 extern PetscClassId IGA_CLASSID;
 
 extern PetscErrorCode IGASetUp(IGA iga);
 
+extern PetscErrorCode IGACreateDM(IGA iga,PetscInt dof,DM *dm);
+extern PetscErrorCode IGACreateDofDM(IGA iga,DM *dm_dof);
+extern PetscErrorCode IGACreateGeomDM(IGA iga,DM *dm_geom);
+
 extern PetscErrorCode IGAGetComm(IGA iga,MPI_Comm *comm);
 extern PetscErrorCode IGAGetDofDM(IGA iga,DM *dm_dof);
 extern PetscErrorCode IGAGetGeomDM(IGA iga,DM *dm_geom);
 
   PetscInt start[3];
   PetscInt width[3];
-  PetscInt *mapping;   /*   [nen] */
-
-  PetscInt     nfix;
-  PetscInt    *ifix;
-  PetscScalar *vfix;
-  PetscScalar *xfix;
+  PetscInt  *mapping;  /*   [nen]        */
+  PetscReal *geometry; /*   [nen][dim+1] */
 
   PetscReal *point;    /*   [nqp][dim]                */
   PetscReal *weight;   /*   [nqp]                     */
   IGA      parent;
   IGAPoint iterator;
 
+  PetscInt     nfix;
+  PetscInt    *ifix;
+  PetscScalar *vfix;
+  PetscScalar *xfix;
+
   PetscInt    nvec;
   PetscScalar *wvec[8];
   PetscInt    nmat;
 
 extern PetscErrorCode IGAElementBegin(IGAElement element);
 extern PetscBool      IGAElementNext(IGAElement element);
+extern PetscErrorCode IGAElementEnd(IGAElement element);
 
 extern PetscErrorCode IGAElementBuildFix(IGAElement element);
 extern PetscErrorCode IGAElementBuildMap(IGAElement element);
   }
   ierr = IGAElementCreate(&iga->iterator);CHKERRQ(ierr);
 
+  iga->geometry = PETSC_NULL;
   iga->rational = PETSC_FALSE;
-  iga->geometry = PETSC_FALSE;
-  iga->dm_geom  = 0;
-  iga->dm_dof   = 0;
+  iga->vec_geom = PETSC_NULL;
+  iga->dm_geom  = PETSC_NULL;
+  iga->dm_dof   = PETSC_NULL;
 
   PetscFunctionReturn(0);
 }
   }
   ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
 
+  ierr = VecDestroy(&iga->vec_geom);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->dm_geom);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->dm_dof);CHKERRQ(ierr);
 
   {
     MPI_Comm comm;
     PetscInt i,dim,dof;
+    PetscBool geometry = iga->vec_geom ? PETSC_TRUE : PETSC_FALSE;
+    PetscBool rational = geometry ? iga->rational : PETSC_FALSE;
     ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
     ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
     ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
 
-    ierr = PetscViewerASCIIPrintf(viewer,"IGA: dimension=%D  dofs/node=%D\n",dim,dof);CHKERRQ(ierr);
+    ierr = PetscViewerASCIIPrintf(viewer,"IGA: dimension=%D  dofs/node=%D  geometry=%s  rational=%s\n",
+                                  dim,dof,geometry?"yes":"no",rational?"yes":"no");CHKERRQ(ierr);
     for (i=0; i<dim; i++) {
       IGAAxis *AX = iga->axis;
       IGARule *QR = iga->rule;
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGACreateDM"
-static PetscErrorCode IGACreateDM(IGA iga,PetscInt dof,DM *_dm)
+PetscErrorCode IGACreateDM(IGA iga,PetscInt dof,DM *_dm)
 {
   PetscInt         i,dim;
   PetscInt         procs[3]   = {-1,-1,-1};
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGACreateDofDM"
-static PetscErrorCode IGACreateDofDM(IGA iga,DM *dm_dof)
+PetscErrorCode IGACreateDofDM(IGA iga,DM *dm_dof)
 {
   PetscInt         dof;
   PetscErrorCode   ierr;
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGACreateGeomDM"
-static PetscErrorCode IGACreateGeomDM(IGA iga,DM *dm_geom)
+PetscErrorCode IGACreateGeomDM(IGA iga,DM *dm_geom)
 {
   PetscInt       dim;
   PetscErrorCode ierr;
 end subroutine IGA_Quadrature_1D
 
 subroutine IGA_ShapeFuns_1D(&
-     rational,geometry,  &
+     geometry,rational,  &
      inq,ina,ind,iJ,iN,  &
      Cw,                 &
      detJ,J,N0,N1,N2,N3) &
   use ISO_C_BINDING, only: C_FLOAT, C_DOUBLE
   implicit none
   integer(kind=C_INT   ), parameter        :: dim = 1
+  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: rational
-  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: inq, ina, ind
   real   (kind=C_DOUBLE), intent(in)  :: iJ, iN(0:ind,ina,inq)
   real   (kind=C_DOUBLE), intent(in)  :: Cw(dim+1,ina)
 
   real(kind=C_DOUBLE) :: C(dim,ina)
   real(kind=C_DOUBLE) :: w(    ina)
+  if (geometry /= 0) then
+     C = Cw(1:dim,:)
+  end if
   if (rational /= 0) then
      w = Cw(dim+1,:) 
   end if
-  if (geometry /= 0) then
-     C = Cw(1:dim,:)
-  end if
 
   nd = max(1,min(ind,3))
   na = ina
 end subroutine IGA_Quadrature_2D
 
 subroutine IGA_ShapeFuns_2D(&
-     rational,geometry,  &
+     geometry,rational,  &
      inq,ina,ind,iJ,iN,  &
      jnq,jna,jnd,jJ,jN,  &
      Cw,                 &
   use ISO_C_BINDING, only: C_FLOAT, C_DOUBLE
   implicit none
   integer(kind=C_INT   ), parameter        :: dim = 2
+  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: rational
-  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: inq, ina, ind
   integer(kind=C_INT   ), intent(in),value :: jnq, jna, jnd
   real   (kind=C_DOUBLE), intent(in)  :: iJ, iN(0:ind,ina,inq)
 
   real(kind=C_DOUBLE) :: C(dim,ina,jna)
   real(kind=C_DOUBLE) :: w(    ina,jna)
+  if (geometry /= 0) then
+     C = Cw(1:dim,:,:)
+  end if
   if (rational /= 0) then
      w = Cw(dim+1,:,:) 
   end if
-  if (geometry /= 0) then
-     C = Cw(1:dim,:,:)
-  end if
 
   nd = max(1,min(ind,jnd,3))
   na = ina*jna
 end subroutine IGA_Quadrature_3D
 
 subroutine IGA_ShapeFuns_3D(&
-     rational,geometry,  &
+     geometry,rational,  &
      inq,ina,ind,iJ,iN,  &
      jnq,jna,jnd,jJ,jN,  &
      knq,kna,knd,kJ,kN,  &
   use ISO_C_BINDING, only: C_FLOAT, C_DOUBLE
   implicit none
   integer(kind=C_INT   ), parameter        :: dim = 3
+  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: rational
-  integer(kind=C_INT   ), intent(in),value :: geometry
   integer(kind=C_INT   ), intent(in),value :: inq, ina, ind
   integer(kind=C_INT   ), intent(in),value :: jnq, jna, jnd
   integer(kind=C_INT   ), intent(in),value :: knq, kna, knd
 
   real(kind=C_DOUBLE) :: C(dim,ina,jna,kna)
   real(kind=C_DOUBLE) :: w(    ina,jna,kna)
+  if (geometry /= 0) then
+     C = Cw(1:dim,:,:,:)
+  end if
   if (rational /= 0) then
      w = Cw(dim+1,:,:,:) 
   end if
-  if (geometry /= 0) then
-     C = Cw(1:dim,:,:,:)
-  end if
 
   nd = max(1,min(ind,jnd,knd,3))
   na = ina*jna*kna
   if (--element->refct > 0) PetscFunctionReturn(0);
 
   ierr = PetscFree(element->mapping);CHKERRQ(ierr);
-
-  ierr = PetscFree(element->ifix);CHKERRQ(ierr);
-  ierr = PetscFree(element->vfix);CHKERRQ(ierr);
-  ierr = PetscFree(element->xfix);CHKERRQ(ierr);
+  ierr = PetscFree(element->geometry);CHKERRQ(ierr);
 
   ierr = PetscFree(element->point);CHKERRQ(ierr);
   ierr = PetscFree(element->weight);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[3]);CHKERRQ(ierr);
   ierr = IGAElementFreeWork(element);CHKERRQ(ierr);
   ierr = IGAPointDestroy(&element->iterator);CHKERRQ(ierr);
+
+  ierr = PetscFree(element->ifix);CHKERRQ(ierr);
+  ierr = PetscFree(element->vfix);CHKERRQ(ierr);
+  ierr = PetscFree(element->xfix);CHKERRQ(ierr);
+
   ierr = PetscFree(element);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
     PetscInt nen = element->nen;
 
     ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*(dim+1),PetscReal,&element->geometry);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim,PetscReal,&element->point);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp,PetscReal,&element->weight);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp,PetscReal,&element->detJ);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
   }
+  element->iterator->parent = element;
+  ierr = IGAPointSetUp(element->iterator);CHKERRQ(ierr);
   { /* */
     PetscInt nen = element->nen;
     PetscInt dof = element->dof;
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->vfix);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->xfix);CHKERRQ(ierr);
   }
-
-  element->iterator->parent = element;
-  ierr = IGAPointSetUp(element->iterator);CHKERRQ(ierr);
-
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "IGAElementBegin"
 PetscErrorCode IGAElementBegin(IGAElement element)
 {
-  IGA iga;
+  IGA            iga;
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   iga = element->parent;
   if (PetscUnlikely(!iga->setup))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetUp() first");
+  if (iga->vec_geom) {
+    ierr = VecGetArrayRead(iga->vec_geom,&iga->geometry);CHKERRQ(ierr);
+  }
   element->index = -1;
   PetscFunctionReturn(0);
 }
     index = (index - coord) / size;
     ID[i] = coord + range[0];
   }
+  IGAElementBuildMap(element);
   IGAElementBuildFix(element);
-  IGAElementBuildMap(element);
   return PETSC_TRUE;
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAElementEnd"
+PetscErrorCode IGAElementEnd(IGAElement element)
+{
+  IGA            iga;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  iga = element->parent;
+  if (iga->vec_geom) {
+    ierr = VecRestoreArrayRead(iga->vec_geom,&iga->geometry);CHKERRQ(ierr);
+  }
+  element->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAElementGetParent"
 PetscErrorCode IGAElementGetParent(IGAElement element,IGA *parent)
 {
   PetscValidPointer(element,1);
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  {
+  { /* */
     IGABasis *BD = element->parent->basis;
     PetscInt *ID = element->ID;
     PetscInt *start = element->start;
       }
     }
   }
+  if (element->parent->geometry) {
+    const PetscScalar *arrayX = element->parent->geometry;
+    PetscReal *Cw = element->geometry;
+    PetscInt *mapping = element->mapping;
+    PetscInt nen = element->nen;
+    PetscInt dim = element->dim;
+    PetscInt a,i,pos=0;
+    for (a=0; a<nen; a++) {
+      const PetscScalar *X = arrayX + mapping[a]*(dim+1);
+      for (i=0; i<(dim+1); i++)
+        Cw[pos++] = PetscRealPart(X[i]);
+    }
+  }
   PetscFunctionReturn(0);
 }
 
 EXTERN_C_BEGIN
 extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
-                             const PetscScalar[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+                             const PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
-                             const PetscScalar[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+                             const PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
                              PetscInt,PetscInt,PetscInt,const PetscReal[],const PetscReal[],
-                             const PetscScalar[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
+                             const PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 EXTERN_C_END
 
 #define IGA_ShapeFuns_ARGS(ID,BD,i) \
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   {
-    PetscBool rational = element->parent->rational;
-    PetscBool geometry = element->parent->geometry;
+    PetscBool geometry = element->parent->geometry ? PETSC_TRUE : PETSC_FALSE;
+    PetscBool rational = geometry ? element->parent->rational : PETSC_FALSE;
     IGABasis *BD  = element->parent->basis;
     PetscInt *ID  = element->ID;
-    PetscScalar *Cw = NULL; /* XXX */
+    PetscReal *Cw = element->geometry;
     PetscReal *dJ = element->detJ;
     PetscReal *J  = element->jacobian;
     PetscReal **N = element->shape;
     switch (element->dim) {
-    case 3: IGA_ShapeFuns_3D(rational,geometry,
+    case 3: IGA_ShapeFuns_3D(geometry,rational,
                              IGA_ShapeFuns_ARGS(ID,BD,0),
                              IGA_ShapeFuns_ARGS(ID,BD,1),
                              IGA_ShapeFuns_ARGS(ID,BD,2),
                              Cw,dJ,J,N[0],N[1],N[2],N[3]); break;
-    case 2: IGA_ShapeFuns_2D(rational,geometry,
+    case 2: IGA_ShapeFuns_2D(geometry,rational,
                              IGA_ShapeFuns_ARGS(ID,BD,0),
                              IGA_ShapeFuns_ARGS(ID,BD,1),
                              Cw,dJ,J,N[0],N[1],N[2],N[3]); break;
-    case 1: IGA_ShapeFuns_1D(rational,geometry,
+    case 1: IGA_ShapeFuns_1D(geometry,rational,
                              IGA_ShapeFuns_ARGS(ID,BD,0),
                              Cw,dJ,J,N[0],N[1],N[2],N[3]); break;
     }
   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
 
   { /* */
-    PetscInt i,dim;
+    PetscInt i,buf[3];
+    PetscInt kind,dim;
     if (!skipheader) {
       PetscInt classid;
       ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr);
       if (classid != IGA_FILE_CLASSID) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Not an IGA in file");
     }
-    ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryRead(viewer,&kind,1,PETSC_INT);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryRead(viewer,&dim, 1,PETSC_INT);CHKERRQ(ierr);
     ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
     for (i=0; i<dim; i++) {
       IGAAxis   axis;
       PetscBool periodic;
-      PetscInt  p,n,m,buf[3];
+      PetscInt  p,m;
       PetscReal *U;
-      ierr = PetscViewerBinaryRead(viewer,&buf,3,PETSC_INT);CHKERRQ(ierr);
+      ierr = PetscViewerBinaryRead(viewer,buf,3,PETSC_INT);CHKERRQ(ierr);
       periodic = buf[0] ? PETSC_TRUE : PETSC_FALSE;
-      p = buf[1]; n = buf[2]; m = n+p+1;
+      p = buf[1];
+      m = buf[2]-1;
       ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
       ierr = IGAAxisSetPeriodic(axis,periodic);CHKERRQ(ierr);
       ierr = IGAAxisSetOrder(axis,p);CHKERRQ(ierr);CHKERRQ(ierr);
       ierr = IGAAxisGetKnots(axis,0,&U);CHKERRQ(ierr);CHKERRQ(ierr);
       ierr = PetscViewerBinaryRead(viewer,U,m+1,PETSC_REAL);CHKERRQ(ierr);
     }
+    if (kind) {
+      DM  dm_geom;
+      Vec vec_geom_global,vec_geom_local;
+      PetscScalar min_w,max_w;
+      ierr = IGACreateGeomDM(iga,&dm_geom);CHKERRQ(ierr);
+      ierr = DMCreateGlobalVector(dm_geom,&vec_geom_global);CHKERRQ(ierr);
+      ierr = VecLoad(vec_geom_global,viewer);CHKERRQ(ierr);
+      ierr = DMCreateLocalVector(dm_geom,&vec_geom_local);CHKERRQ(ierr);
+      ierr = DMGlobalToLocalBegin(dm_geom,vec_geom_global,INSERT_VALUES,vec_geom_local);CHKERRQ(ierr);
+      ierr = DMGlobalToLocalEnd  (dm_geom,vec_geom_global,INSERT_VALUES,vec_geom_local);CHKERRQ(ierr);
+
+      ierr = VecDuplicate(vec_geom_local,&iga->vec_geom);CHKERRQ(ierr);
+      ierr = VecCopy(vec_geom_local,iga->vec_geom);CHKERRQ(ierr);
+      ierr = VecStrideMin(vec_geom_global,dim,PETSC_NULL,&min_w);CHKERRQ(ierr);
+      ierr = VecStrideMax(vec_geom_global,dim,PETSC_NULL,&max_w);CHKERRQ(ierr);
+      iga->rational = (PetscAbs(max_w-min_w) > 100*PETSC_MACHINE_EPSILON) ? PETSC_TRUE : PETSC_FALSE;
+
+      ierr = VecDestroy(&vec_geom_global);CHKERRQ(ierr);
+      ierr = VecDestroy(&vec_geom_local);CHKERRQ(ierr);
+      ierr = DMDestroy(&dm_geom);CHKERRQ(ierr);
+    }
   }
   PetscFunctionReturn(0);
 }
   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
 
   { /* */
-    PetscInt i,dim;
+    PetscInt i=0,buf[3];
+    PetscInt kind,dim;
+    kind = iga->vec_geom ? (1 + (PetscInt)iga->rational) : 0;
     ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-    if (!skipheader) {
-      PetscInt classid = IGA_FILE_CLASSID;
-      ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
-    }
-    ierr = PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr);
+    if (!skipheader) buf[i++] = IGA_FILE_CLASSID;
+    buf[i++] = kind; buf[i++] = dim;
+    ierr = PetscViewerBinaryWrite(viewer,buf,i,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
     for (i=0; i<dim; i++) {
       IGAAxis   axis;
       PetscBool periodic;
-      PetscInt  p,n,m,buf[3];
+      PetscInt  p,m;
       PetscReal *U;
       ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
       ierr = IGAAxisGetPeriodic(axis,&periodic);CHKERRQ(ierr);
       ierr = IGAAxisGetKnots(axis,&m,&U);CHKERRQ(ierr);
       buf[0] = periodic;
       buf[1] = p;
-      buf[2] = n = m-p-1;
+      buf[2] = m+1;
       ierr = PetscViewerBinaryWrite(viewer,buf,3,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
       ierr = PetscViewerBinaryWrite(viewer,U,m+1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
     }
+    if (iga->vec_geom) {
+      DM  dm_geom;
+      Vec vec_geom_global,vec_geom_local;
+      ierr = IGACreateGeomDM(iga,&dm_geom);CHKERRQ(ierr);
+      ierr = DMCreateLocalVector(dm_geom,&vec_geom_local);CHKERRQ(ierr);
+      ierr = VecCopy(iga->vec_geom,vec_geom_local);CHKERRQ(ierr);
+      ierr = DMCreateGlobalVector(dm_geom,&vec_geom_global);CHKERRQ(ierr);
+      ierr = DMLocalToGlobalBegin(dm_geom,vec_geom_local,INSERT_VALUES,vec_geom_global);CHKERRQ(ierr);
+      ierr = DMLocalToGlobalEnd  (dm_geom,vec_geom_local,INSERT_VALUES,vec_geom_global);CHKERRQ(ierr);
+      ierr = VecView(vec_geom_global,viewer);CHKERRQ(ierr);
+      ierr = VecDestroy(&vec_geom_global);CHKERRQ(ierr);
+      ierr = VecDestroy(&vec_geom_local);CHKERRQ(ierr);
+      ierr = DMDestroy(&dm_geom);CHKERRQ(ierr);
+    }
     ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidCharPointer(filename,2);
   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
-  ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
+  ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr);
+  ierr = PetscViewerSetType(viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
+  ierr = PetscViewerBinarySkipInfo(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerFileSetMode(viewer,FILE_MODE_READ);CHKERRQ(ierr);
+  ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr);
   ierr = IGALoad(iga,viewer);CHKERRQ(ierr);
   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
   PetscFunctionReturn(0);
     ierr = IGAElementAssembleMat(element,A,matA);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
   }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
 
   ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
   ierr = MatAssemblyEnd  (matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

src/petigapoint.c

   element = point->parent;
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  point->index = -1;
   ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
   ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
   PetscFunctionReturn(0);
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
 
   /* Restore local vector U */
   ierr = VecRestoreArrayRead(localU,&arrayU);
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);;CHKERRQ(ierr);
   }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
 
   /* Restore local vector U */
   ierr = VecRestoreArrayRead(localU,&arrayU);
     ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
     ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
   }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
 
   /* Restore array to the local vectors */
   ierr = VecRestoreArrayRead(localV,&arrayV);CHKERRQ(ierr);
     ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
     ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
   }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
 
   /* Restore array to the vectors */
   ierr = VecRestoreArrayRead(localU,&arrayU);CHKERRQ(ierr);