Commits

Lisandro Dalcin committed a0fe97f

Add support for additional field data, aka property

Comments (0)

Files changed (6)

demo/ClassicalShell.c

   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
-  ierr = IGASetSpatialDim(iga,3);CHKERRQ(ierr);
+  ierr = IGASetGeometryDim(iga,3);CHKERRQ(ierr);
   ierr = IGASetDof(iga,5);CHKERRQ(ierr); // dofs = {ux,uy,uz,psix,psiy}
   ierr = IGARead(iga,filename);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   char       **fieldname;
 
   PetscInt  dim;   /* parametric dimension of the function space*/
-  PetscInt  nsd;   /* spatial dimension of the geometry */
   PetscInt  dof;   /* number of degrees of freedom per node */
   PetscInt  order; /* maximum derivative order */
 
   IGABasis    node_basis[3];
   IGAElement  node_iterator;
 
+
+  PetscInt    geometry;
+  PetscBool   rational;
+  PetscInt    property;
+  PetscReal   *geometryX;
+  PetscReal   *rationalW;
+  PetscScalar *propertyA;
+
   PetscInt  proc_sizes[3];
   PetscInt  proc_ranks[3];
 
   Vec       elem_vec;
   DM        elem_dm;
 
-  PetscBool geometry;
-  PetscBool rational;
-  PetscReal *geometryX;
-  PetscReal *geometryW;
   PetscInt  geom_sizes[3];
   PetscInt  geom_lstart[3];
   PetscInt  geom_lwidth[3];
   PetscInt  geom_gstart[3];
   PetscInt  geom_gwidth[3];
-  Vec       geom_vec;
   DM        geom_dm;
 
   PetscInt  node_sizes[3];
 PETSC_EXTERN PetscErrorCode IGARead(IGA iga,const char filename[]);
 PETSC_EXTERN PetscErrorCode IGAWrite(IGA iga,const char filename[]);
 
+PETSC_EXTERN PetscErrorCode IGASetGeometryDim(IGA iga,PetscInt dim);
+PETSC_EXTERN PetscErrorCode IGAGetGeometryDim(IGA iga,PetscInt *dim);
 PETSC_EXTERN PetscErrorCode IGALoadGeometry(IGA iga,PetscViewer viewer);
 PETSC_EXTERN PetscErrorCode IGASaveGeometry(IGA iga,PetscViewer viewer);
 
+PETSC_EXTERN PetscErrorCode IGASetPropertyDim(IGA iga,PetscInt dim);
+PETSC_EXTERN PetscErrorCode IGAGetPropertyDim(IGA iga,PetscInt *dim);
+PETSC_EXTERN PetscErrorCode IGALoadProperty(IGA iga,PetscViewer viewer);
+PETSC_EXTERN PetscErrorCode IGASaveProperty(IGA iga,PetscViewer viewer);
+
 PETSC_EXTERN PetscErrorCode IGALoadVec(IGA iga,Vec vec,PetscViewer viewer);
 PETSC_EXTERN PetscErrorCode IGASaveVec(IGA iga,Vec vec,PetscViewer viewer);
 PETSC_EXTERN PetscErrorCode IGAReadVec(IGA iga,Vec vec,const char filename[]);
 
 PETSC_EXTERN PetscErrorCode IGASetDim(IGA iga,PetscInt dim);
 PETSC_EXTERN PetscErrorCode IGAGetDim(IGA iga,PetscInt *dim);
-PETSC_EXTERN PetscErrorCode IGASetSpatialDim(IGA iga,PetscInt nsd);
-PETSC_EXTERN PetscErrorCode IGAGetSpatialDim(IGA iga,PetscInt *nsd);
 PETSC_EXTERN PetscErrorCode IGASetDof(IGA iga,PetscInt dof);
 PETSC_EXTERN PetscErrorCode IGAGetDof(IGA iga,PetscInt *dof);
 PETSC_EXTERN PetscErrorCode IGASetFieldName(IGA iga,PetscInt field,const char name[]);
   PetscInt dof;
   PetscInt dim;
   PetscInt nsd;
+  PetscInt npd;
+
   IGABasis *BD;
 
-  PetscInt  *mapping;  /*   [nen]      */
-
-  PetscBool geometry;
-  PetscBool rational;
-  PetscReal *geometryX;/*   [nen][nsd] */
-  PetscReal *geometryW;/*   [nen]      */
+  PetscInt    *mapping;   /*[nen]      */
+  PetscBool   geometry;
+  PetscReal   *geometryX; /*[nen][nsd] */
+  PetscBool   rational;
+  PetscReal   *rationalW; /*[nen]      */
+  PetscBool   property;
+  PetscScalar *propertyA; /*[nen][npd] */
 
   PetscReal *weight;   /*   [nqp]                     */
   PetscReal *detJac;   /*   [nqp]                     */
 
 PETSC_EXTERN PetscErrorCode IGAElementBuildMapping(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildGeometry(IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAElementBuildProperty(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildQuadrature(IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAElementBuildShapeFuns(IGAElement element);
 
   PetscInt dof;
   PetscInt dim;
   PetscInt nsd;
+  PetscInt npd;
 
-  PetscReal *weight;   /*      */
-  PetscReal *detJac;   /*      */
+  PetscReal   *geometry;/*  [nen][nsd] */
+  PetscScalar *property;/*  [nen][npd] */
 
+  PetscReal *weight;   /*   [1]   */
+  PetscReal *detJac;   /*   [1]   */
   PetscReal *point;    /*   [dim] */
   PetscReal *scale;    /*   [dim] */
   PetscReal *basis[4]; /*0: [nen] */
                        /*2: [nen][dim][dim] */
                        /*3: [nen][dim][dim][dim] */
 
-  PetscReal *geometry; /*   [nen][nsd] */
   PetscReal *detX;     /*   [1] */
   PetscReal *gradX[2]; /*0: [nsd][dim] */
                        /*1: [dim][nsd] */
   ierr = VecDestroy(&iga->elem_vec);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->elem_dm);CHKERRQ(ierr);
   /* geometry */
-  iga->geometry = PETSC_FALSE;
+  iga->geometry = 0;
   iga->rational = PETSC_FALSE;
+  iga->property = 0;
   ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
-  ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
-  ierr = VecDestroy(&iga->geom_vec);CHKERRQ(ierr);
+  ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
+  ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->geom_dm);CHKERRQ(ierr);
   /* node */
   ierr = AODestroy(&iga->ao);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGASetSpatialDim"
-/*@
-   IGASetSpatialDim - Sets the dimension of the geometry
-
-   Logically Collective on IGA
-
-   Input Parameters:
-+  iga - the IGA context
--  dim - the dimension of the geometry
-
-   Level: normal
-
-.keywords: IGA, dimension
-@*/
-PetscErrorCode IGASetSpatialDim(IGA iga,PetscInt nsd)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidLogicalCollectiveInt(iga,nsd,2);
-  if (nsd < 1 || nsd > 3)
-    SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
-             "Number of space dimensions must be in range [1,3], got %D",nsd);
-  if (iga->nsd > 0 && iga->nsd != nsd)
-    SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
-             "Cannot change IGA dim from %D after it was set to %D",iga->nsd,nsd);
-  iga->nsd = nsd;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGAGetSpatialDim"
-PetscErrorCode IGAGetSpatialDim(IGA iga,PetscInt *nsd)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidPointer(nsd,2);
-  *nsd = iga->nsd;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
 #define __FUNCT__ "IGASetDof"
 /*@
    IGASetDof - Sets the number of degrees of freedom per basis
   ierr = VecDestroy(&iga->elem_vec);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->elem_dm);CHKERRQ(ierr);
   /* geometry */
-  iga->geometry = PETSC_FALSE;
+  iga->geometry = 0;
   iga->rational = PETSC_FALSE;
+  iga->property = 0;
   ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
-  ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
-  ierr = VecDestroy(&iga->geom_vec);CHKERRQ(ierr);
+  ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
+  ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->geom_dm);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 
   ierr = PetscFree(element->mapping);CHKERRQ(ierr);
   ierr = PetscFree(element->geometryX);CHKERRQ(ierr);
-  ierr = PetscFree(element->geometryW);CHKERRQ(ierr);
+  ierr = PetscFree(element->rationalW);CHKERRQ(ierr);
+  ierr = PetscFree(element->propertyA);CHKERRQ(ierr);
 
   ierr = PetscFree(element->weight);CHKERRQ(ierr);
   ierr = PetscFree(element->detJac);CHKERRQ(ierr);
 
   element->dof = iga->dof;
   element->dim = iga->dim;
-  element->nsd = iga->nsd ? iga->nsd : iga->dim;
+  element->nsd = iga->geometry ? iga->geometry : iga->dim;
+  element->npd = iga->property ? iga->property : 0;
   if (!element->collocation) {
     start = iga->elem_start;
     width = iga->elem_width;
   { /* */
     PetscInt dim = element->dim;
     PetscInt nsd = element->nsd;
+    PetscInt npd = element->npd;
     PetscInt nen = element->nen;
     PetscInt nqp = element->nqp;
 
     ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*nsd,PetscReal,&element->geometryX);CHKERRQ(ierr);
-    ierr = PetscMalloc1(nen    ,PetscReal,&element->geometryW);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen,PetscReal,&element->rationalW);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*npd,PetscScalar,&element->propertyA);CHKERRQ(ierr);
 
     ierr = PetscMalloc1(nqp,PetscReal,&element->weight);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp,PetscReal,&element->detJac);CHKERRQ(ierr);
 #define CHKERRRETURN(n,r) do{if(PetscUnlikely(n)){CHKERRCONTINUE(n);return(r);}}while(0)
   ierr = IGAElementBuildMapping(element);  CHKERRRETURN(ierr,PETSC_FALSE);
   ierr = IGAElementBuildGeometry(element); CHKERRRETURN(ierr,PETSC_FALSE);
+  ierr = IGAElementBuildProperty(element); CHKERRRETURN(ierr,PETSC_FALSE);
   ierr = IGAElementBuildFix(element);      CHKERRRETURN(ierr,PETSC_FALSE);
 #undef  CHKERRRETURN
   return PETSC_TRUE;
   (*point)->dof = element->dof;
   (*point)->dim = element->dim;
   (*point)->nsd = element->nsd;
+  (*point)->npd = element->npd;
 
   PetscFunctionReturn(0);
 }
 
  start:
 
+  point->geometry = element->geometryX;
+  point->property = element->propertyA;
+  if (!element->geometry)
+    point->geometry = PETSC_NULL;
+  if (!element->property)
+    point->property = PETSC_NULL;
+  
   point->weight   = element->weight;
   point->detJac   = element->detJac;
 
   point->basis[2] = element->basis[2];
   point->basis[3] = element->basis[3];
 
-  if (element->geometry)
-    point->geometry = element->geometryX;
-  else
-    point->geometry = PETSC_NULL;
   if (element->geometry && dim == nsd) { /* XXX */
     point->detX     = element->detX;
     point->gradX[0] = element->gradX[0];
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
   iga = element->parent;
-  element->geometry = iga->geometry;
-  element->rational = iga->rational;
+  element->geometry = iga->geometry ? PETSC_TRUE : PETSC_FALSE;
+  element->rational = iga->rational ? PETSC_TRUE : PETSC_FALSE;
   if (element->geometry || element->rational) {
     const PetscInt  *map = element->mapping;
     const PetscReal *arrayX = iga->geometryX;
-    const PetscReal *arrayW = iga->geometryW;
+    const PetscReal *arrayW = iga->rationalW;
     PetscReal *X = element->geometryX;
-    PetscReal *W = element->geometryW;
+    PetscReal *W = element->rationalW;
     PetscInt a,nen = element->nen;
     PetscInt i,nsd = element->nsd;
     if (element->geometry)
   PetscFunctionReturn(0);
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementBuildProperty"
+PetscErrorCode IGAElementBuildProperty(IGAElement element)
+{
+  IGA iga;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  iga = element->parent;
+  element->property = iga->property ? PETSC_TRUE: PETSC_FALSE;
+  if (element->property) {
+    const PetscInt *map = element->mapping;
+    const PetscScalar *arrayA = iga->propertyA;
+    PetscScalar *A = element->propertyA;
+    PetscInt a,nen = element->nen;
+    PetscInt i,npd = element->npd;
+    for (a=0; a<nen; a++)
+      for (i=0; i<npd; i++)
+        A[i + a*npd] = arrayA[i + map[a]*npd];
+  }
+  PetscFunctionReturn(0);
+}
+
 EXTERN_C_BEGIN
 extern void IGA_Quadrature_1D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                               PetscReal[],PetscReal[],PetscReal*,PetscReal[]);
     PetscReal **N = element->basis;
     switch (element->dim) {
     case 3: IGA_BasisFuns_3D(order,element->rational,
-                             element->geometryW,
+                             element->rationalW,
                              IGA_BasisFuns_ARGS(ID,BD,0),
                              IGA_BasisFuns_ARGS(ID,BD,1),
                              IGA_BasisFuns_ARGS(ID,BD,2),
                              N[0],N[1],N[2],N[3]); break;
     case 2: IGA_BasisFuns_2D(order,element->rational,
-                             element->geometryW,
+                             element->rationalW,
                              IGA_BasisFuns_ARGS(ID,BD,0),
                              IGA_BasisFuns_ARGS(ID,BD,1),
                              N[0],N[1],N[2],N[3]); break;
     case 1: IGA_BasisFuns_1D(order,element->rational,
-                             element->geometryW,
+                             element->rationalW,
                              IGA_BasisFuns_ARGS(ID,BD,0),
                              N[0],N[1],N[2],N[3]); break;
     }
     switch (dim) {
     case 2: IGA_BoundaryArea_2D(shape,dir,side,
                                 element->geometry,element->geometryX,
-                                element->rational,element->geometryW,
+                                element->rational,element->rationalW,
                                 nqp[0],W[0],nen[0],ndr[0],N[0],
                                 &dS); break;
     case 3: IGA_BoundaryArea_3D(shape,dir,side,
                                 element->geometry,element->geometryX,
-                                element->rational,element->geometryW,
+                                element->rational,element->rationalW,
                                 nqp[0],W[0],nen[0],ndr[0],N[0],
                                 nqp[1],W[1],nen[1],ndr[1],N[1],
                                 &dS);break;
      integer(kind=IGA_INTEGER_KIND) :: dof
      integer(kind=IGA_INTEGER_KIND) :: dim
      integer(kind=IGA_INTEGER_KIND) :: nsd
+     integer(kind=IGA_INTEGER_KIND) :: npd
+
+     type(C_PTR) :: geometry
+     type(C_PTR) :: property
 
      type(C_PTR) :: weight
      type(C_PTR) :: detJac
      type(C_PTR) :: point
      type(C_PTR) :: scale
      type(C_PTR) :: basis(0:3)
-     type(C_PTR) :: geometry
      type(C_PTR) :: detX
      type(C_PTR) :: gradX(0:1)
      type(C_PTR) :: shape(0:3)
 
   contains
 
+    pure function IGA_Geometry(p) result(X)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: X(:)
+      call c2f(p%geometry,X,(/p%nsd,p%nen/))
+    end function IGA_Geometry
+
+    pure function IGA_Property(p) result(A)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: A(:)
+      call c2f(p%property,A,(/p%npd,p%nen/))
+    end function IGA_Property
+
     function IGA_GeomMap(p) result (X)
       implicit none
       type(IGAPoint), intent(in) :: p
 {
   PetscBool      isbinary;
   PetscBool      skipheader;
-  PetscInt       descr;
+  PetscBool      geometry;
+  PetscBool      property;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   if (!isbinary) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Only for binary viewers");
   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
 
-  ierr = IGAReset(iga);CHKERRQ(ierr);
-
   if (!skipheader) {
     PetscInt classid = 0;
     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,&descr,1,PETSC_INT);CHKERRQ(ierr);
-  if (descr >= 0) { /* */
+  { /* */
+    PetscInt info = 0;
+    ierr = PetscViewerBinaryRead(viewer,&info,1,PETSC_INT);CHKERRQ(ierr);
+    geometry = (info & 0x1) ? PETSC_TRUE : PETSC_FALSE; 
+    property = (info & 0x2) ? PETSC_TRUE : PETSC_FALSE;
+  }
+  ierr = IGAReset(iga);CHKERRQ(ierr);
+  { /* */
     PetscInt i,dim;
     ierr = PetscViewerBinaryRead(viewer,&dim, 1,PETSC_INT);CHKERRQ(ierr);
     ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
       ierr = IGAAxisGetKnots(axis,PETSC_NULL,&U);CHKERRQ(ierr);CHKERRQ(ierr);
       ierr = PetscViewerBinaryRead(viewer,U,m,PETSC_REAL);CHKERRQ(ierr);
     }
-    ierr = IGASetUp_Basic(iga);CHKERRQ(ierr);
   }
-  if (PetscAbs(descr) >= 1) { /* */
-    PetscInt nsd;
-    ierr = PetscViewerBinaryRead(viewer,&nsd,1,PETSC_INT);CHKERRQ(ierr);
-    ierr = IGASetSpatialDim(iga,nsd);CHKERRQ(ierr);
+  ierr = IGASetUp_Basic(iga);CHKERRQ(ierr);
+  if (geometry) { /* */
+    PetscInt dim;
+    ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr);
+    ierr = IGASetGeometryDim(iga,dim);CHKERRQ(ierr);
     ierr = IGALoadGeometry(iga,viewer);CHKERRQ(ierr);
   }
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
-#if 0
-  /* XXX waiting implementation ... */
-  if (PetscAbs(descr) >= 2) { /* */
-    PetscInt npd,ncd;
-    ierr = PetscViewerBinaryRead(viewer,&npd,1,PETSC_INT);CHKERRQ(ierr);
-    ierr = PetscViewerBinaryRead(viewer,&ncd,1,PETSC_INT);CHKERRQ(ierr);
-    ierr = IGASetDataDim(iga,npd,ncd);CHKERRQ(ierr);
-    if (npd > 0) {ierr = IGALoadPointData(iga,viewer);CHKERRQ(ierr);}
-    if (ncd > 0) {ierr = IGALoadCellData(iga,viewer);CHKERRQ(ierr);}
+  if (property) { /* */
+    PetscInt dim;
+    ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr);
+    ierr = IGASetPropertyDim(iga,dim);CHKERRQ(ierr);
+    ierr = IGALoadProperty(iga,viewer);CHKERRQ(ierr);
   }
-#endif
   PetscFunctionReturn(0);
 }
 
     ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
   }
   { /* */
-    PetscInt descr = iga->geometry ? 1 : 0;
-    ierr = PetscViewerBinaryWrite(viewer,&descr,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
+    PetscInt info = 0;
+    if (iga->geometry) info |= 0x1;
+    if (iga->property) info |= 0x2;
+    ierr = PetscViewerBinaryWrite(viewer,&info,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
   }
   { /* */
     PetscInt i,dim;
       ierr = PetscViewerBinaryWrite(viewer,U,m+1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
     }
   }
-  if (iga->geometry) {
-    PetscInt nsd;
-    ierr = IGAGetSpatialDim(iga,&nsd);CHKERRQ(ierr);
-    ierr = PetscViewerBinaryWrite(viewer,&nsd,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
+  if (iga->geometry) { /* */
+    PetscInt dim;
+    ierr = IGAGetGeometryDim(iga,&dim);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
     ierr = IGASaveGeometry(iga,viewer);CHKERRQ(ierr);
   }
+  if (iga->property) { /* */
+    PetscInt dim;
+    ierr = IGAGetPropertyDim(iga,&dim);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = IGASaveProperty(iga,viewer);CHKERRQ(ierr);
+  }
   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGASetGeometryDim"
+/*@
+   IGASetGeometryDim - Sets the dimension of the geometry
+
+   Logically Collective on IGA
+
+   Input Parameters:
++  iga - the IGA context
+-  dim - the dimension of the geometry
+
+   Level: normal
+
+.keywords: IGA, dimension
+@*/
+PetscErrorCode IGASetGeometryDim(IGA iga,PetscInt dim)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidLogicalCollectiveInt(iga,dim,2);
+  if (dim < 1 || dim > 3)
+    SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+             "Number of space dimensions must be in range [1,3], got %D",dim);
+  if (iga->geometry > 0 && iga->geometry != dim)
+    SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+             "Cannot change IGA spatial dim to %D after it was set to %D",dim,iga->geometry);
+  if (iga->geometry == 0) iga->setup = PETSC_FALSE;
+  iga->geometry = dim;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAGetGeometryDim"
+PetscErrorCode IGAGetGeometryDim(IGA iga,PetscInt *dim)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(dim,2);
+  *dim = iga->geometry;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGALoadGeometry"
 PetscErrorCode IGALoadGeometry(IGA iga,PetscViewer viewer)
 {
   if (!isbinary) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Only for binary viewers");
   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
 
-  ierr = IGAGetSpatialDim(iga,&nsd);CHKERRQ(ierr);
+  ierr = IGAGetGeometryDim(iga,&nsd);CHKERRQ(ierr);
   if (nsd < 1) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
-                       "Must call IGASetSpatialDim() first");
+                       "Must call IGASetGeometryDim() first");
   {
     IGA_Grid grid;
     ierr = IGA_NewGridIO(iga,nsd+1,&grid);CHKERRQ(ierr);
   ierr = VecStrideMin(gvec,nsd,PETSC_NULL,&min_w);CHKERRQ(ierr);
   ierr = VecStrideMax(gvec,nsd,PETSC_NULL,&max_w);CHKERRQ(ierr);
 
-  iga->geometry = PETSC_TRUE;
   iga->rational = ((max_w-min_w)>tol_w) ? PETSC_TRUE : PETSC_FALSE;
   ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
-  ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
+  ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
   {
     PetscInt n,a,i,pos;
     PetscReal *X,*W;
     ierr = VecGetSize(lvec,&n);CHKERRQ(ierr);
     n /= (nsd+1);
     ierr = PetscMalloc1(n*nsd,PetscReal,&iga->geometryX);CHKERRQ(ierr);
-    ierr = PetscMalloc1(n,    PetscReal,&iga->geometryW);CHKERRQ(ierr);
-    X = iga->geometryX; W = iga->geometryW;
+    ierr = PetscMalloc1(n,    PetscReal,&iga->rationalW);CHKERRQ(ierr);
+    X = iga->geometryX; W = iga->rationalW;
     ierr = VecGetArrayRead(lvec,&Xw);CHKERRQ(ierr);
     for (pos=0,a=0; a<n; a++) {
       for (i=0; i<nsd; i++)
   if (!isbinary) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Only for binary viewers");
   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
 
-  ierr = IGAGetSpatialDim(iga,&nsd);CHKERRQ(ierr);
+  ierr = IGAGetGeometryDim(iga,&nsd);CHKERRQ(ierr);
   {
     IGA_Grid grid;
     ierr = IGA_NewGridIO(iga,nsd+1,&grid);CHKERRQ(ierr);
     PetscInt n,a,i,pos;
     PetscScalar *Xw;
     const PetscReal *X = iga->geometryX;
-    const PetscReal *W = iga->geometryW;
+    const PetscReal *W = iga->rationalW;
     ierr = VecGetSize(lvec,&n);CHKERRQ(ierr);
     n /= (nsd+1);
     ierr = VecGetArray(lvec,&Xw);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetPropertyDim"
+/*@
+   IGASetPropertyDim - Sets the dimension of the property
+
+   Logically Collective on IGA
+
+   Input Parameters:
++  iga - the IGA context
+-  dim - the dimension of the property
+
+   Level: normal
+
+.keywords: IGA, dimension
+@*/
+PetscErrorCode IGASetPropertyDim(IGA iga,PetscInt dim)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidLogicalCollectiveInt(iga,dim,2);
+  if (dim < 1)
+    SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+             "Number of properties must be greater than 0, got %D",dim);
+  if (iga->property > 0 && iga->property != dim)
+    SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+             "Cannot change IGA property dim to %D after it was set to %D",dim,iga->property);
+  if (iga->property == 0) iga->setup = PETSC_FALSE;
+  iga->property = dim;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAGetPropertyDim"
+PetscErrorCode IGAGetPropertyDim(IGA iga,PetscInt *dim)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(dim,2);
+  *dim = iga->property;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGALoadProperty"
+PetscErrorCode IGALoadProperty(IGA iga,PetscViewer viewer)
+{
+  PetscBool      isbinary;
+  PetscBool      skipheader;
+  PetscInt       npd;
+  Vec            nvec,gvec,lvec;
+  VecScatter     g2n,g2l;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
+  PetscCheckSameComm(iga,1,viewer,2);
+  if (iga->setupstage < 1) IGACheckSetUp(iga,1);
+
+  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
+  if (!isbinary) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Only for binary viewers");
+  ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
+
+  ierr = IGAGetPropertyDim(iga,&npd);CHKERRQ(ierr);
+  if (npd < 1) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+                       "Must call IGASetPropertyDim() first");
+  {
+    IGA_Grid grid;
+    ierr = IGA_NewGridIO(iga,npd,&grid);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecNatural(grid,iga->vectype,&nvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecGlobal (grid,iga->vectype,&gvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecLocal  (grid,iga->vectype,&lvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetScatterG2N(grid,&g2n);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetScatterG2L(grid,&g2l);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)nvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)gvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)lvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)g2n);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)g2l);CHKERRQ(ierr);
+    ierr = IGA_Grid_Destroy(&grid);CHKERRQ(ierr);
+  }
+  /* viewer -> natural*/
+  if (!skipheader)
+    {ierr = VecLoad(nvec,viewer);CHKERRQ(ierr);}
+  else
+    {ierr = VecLoad_Binary_SkipHeader(nvec,viewer);CHKERRQ(ierr);}
+  /* natural -> global */
+  ierr = VecScatterBegin(g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  ierr = VecScatterEnd  (g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  /* global -> local */
+  ierr = VecScatterBegin(g2l,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd  (g2l,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+
+  ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
+  {
+    PetscInt n; const PetscScalar *A;
+    ierr = VecGetSize(lvec,&n);CHKERRQ(ierr);
+    ierr = PetscMalloc1(n,PetscScalar,&iga->propertyA);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(lvec,&A);CHKERRQ(ierr);
+    ierr = PetscMemcpy(iga->propertyA,A,n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = VecRestoreArrayRead(lvec,&A);CHKERRQ(ierr);
+  }
+
+  ierr = VecScatterDestroy(&g2n);CHKERRQ(ierr);
+  ierr = VecScatterDestroy(&g2l);CHKERRQ(ierr);
+  ierr = VecDestroy(&lvec);CHKERRQ(ierr);
+  ierr = VecDestroy(&gvec);CHKERRQ(ierr);
+  ierr = VecDestroy(&nvec);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASaveProperty"
+PetscErrorCode IGASaveProperty(IGA iga,PetscViewer viewer)
+{
+  PetscBool      isbinary;
+  PetscBool      skipheader;
+  PetscInt       npd;
+  Vec            nvec,gvec,lvec;
+  VecScatter     l2g,g2n;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
+  PetscCheckSameComm(iga,1,viewer,2);
+  if (iga->setupstage < 1) IGACheckSetUp(iga,1);
+  if (!iga->property) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,"No property set");
+
+  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
+  if (!isbinary) SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_ARG_WRONG,"Only for binary viewers");
+  ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
+
+  ierr = IGAGetPropertyDim(iga,&npd);CHKERRQ(ierr);
+  {
+    IGA_Grid grid;
+    ierr = IGA_NewGridIO(iga,npd,&grid);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecNatural(grid,iga->vectype,&nvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecGlobal (grid,iga->vectype,&gvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecLocal  (grid,iga->vectype,&lvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetScatterL2G(grid,&l2g);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetScatterG2N(grid,&g2n);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)nvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)gvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)lvec);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)l2g);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)g2n);CHKERRQ(ierr);
+    ierr = IGA_Grid_Destroy(&grid);CHKERRQ(ierr);
+  }
+  {
+    PetscInt n; PetscScalar *A;
+    ierr = VecGetSize(lvec,&n);CHKERRQ(ierr);
+    ierr = VecGetArray(lvec,&A);CHKERRQ(ierr);
+    ierr = PetscMemcpy(A,iga->propertyA,n*sizeof(PetscScalar));CHKERRQ(ierr);
+    ierr = VecRestoreArray(lvec,&A);CHKERRQ(ierr);
+  }
+  /* local -> global */
+  ierr = VecScatterBegin(l2g,lvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd  (l2g,lvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  /* global -> natural */
+  ierr = VecScatterBegin(g2n,gvec,nvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd  (g2n,gvec,nvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  /* natural -> viewer */
+  ierr = VecView(nvec,viewer);CHKERRQ(ierr);
+
+  ierr = VecScatterDestroy(&g2n);CHKERRQ(ierr);
+  ierr = VecScatterDestroy(&l2g);CHKERRQ(ierr);
+  ierr = VecDestroy(&lvec);CHKERRQ(ierr);
+  ierr = VecDestroy(&gvec);CHKERRQ(ierr);
+  ierr = VecDestroy(&nvec);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGARead"