Commits

Lisandro Dalcin  committed aca0829

Much better implementation of natural<->global scatters

  • Participants
  • Parent commits 3cfad3a

Comments (0)

Files changed (8)

File demo/InputOutput.c

   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
   
   ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
   ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr); /* just for testing */
   ierr = VecDestroy(&vec);CHKERRQ(ierr);
 
   ierr = VecCreate(comm,&vec);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(vec,dof);CHKERRQ(ierr);
   ierr = PetscViewerBinaryOpen(comm,"igavec.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
   ierr = VecLoad(vec,viewer);CHKERRQ(ierr);
   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

File demo/makefile

 runex0a_1:
 	-@${MPIEXEC} -n 1 ./InputOutput ${OPTS} -dim 1 -periodic
 runex0a_2:
-	-@${MPIEXEC} -n 2 ./InputOutput ${OPTS} -dim 2 -periodic
+	-@${MPIEXEC} -n 2 ./InputOutput ${OPTS} -dim 2 -dof 3 -periodic
 runex0a_3:
 	-@${MPIEXEC} -n 3 ./InputOutput ${OPTS} -dim 3 -periodic 1,0,1
 runex0a_4:

File include/petiga.h

   VecScatter g2l,l2g;
   PetscInt   nwork;
   Vec        vwork[16];
-  VecScatter g2n;
   Vec        natural;
+  VecScatter n2g,g2n;
 
 };
 

File src/petiga.c

   ierr = VecScatterDestroy(&iga->g2l);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&iga->l2g);CHKERRQ(ierr);
   ierr = VecDestroy(&iga->natural);CHKERRQ(ierr);
+  ierr = VecScatterDestroy(&iga->n2g);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&iga->g2n);CHKERRQ(ierr);
   while (iga->nwork > 0)
     {ierr = VecDestroy(&iga->vwork[--iga->nwork]);CHKERRQ(ierr);}
       geom_gwidth[i] = 1;
     }
   }
+  /* element */
+  ierr = VecDestroy(&iga->elem_vec);CHKERRQ(ierr);
+  ierr = DMDestroy(&iga->elem_dm);CHKERRQ(ierr);
+  /* geometry */
+  iga->geometry = PETSC_FALSE;
+  iga->rational = PETSC_FALSE;
+  ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
+  ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
+  ierr = VecDestroy(&iga->geom_vec);CHKERRQ(ierr);
+  ierr = DMDestroy(&iga->geom_dm);CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }
 
   ierr = ISLocalToGlobalMappingDestroy(&iga->lgmapb);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&iga->g2l);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&iga->l2g);CHKERRQ(ierr);
+  ierr = VecDestroy(&iga->natural);CHKERRQ(ierr);
+  ierr = VecScatterDestroy(&iga->n2g);CHKERRQ(ierr);
+  ierr = VecScatterDestroy(&iga->g2n);CHKERRQ(ierr);
   while (iga->nwork > 0)
     {ierr = VecDestroy(&iga->vwork[--iga->nwork]);CHKERRQ(ierr);}
   ierr = DMDestroy(&iga->node_dm);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject)iga->lgmap);CHKERRQ(ierr);
     ierr = IGA_Grid_GetLGMapBlock(grid,&iga->lgmapb);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject)iga->lgmapb);CHKERRQ(ierr);
-    /* build global to local and local to global vector scatters */
+    /* build global <-> local vector scatters */
     ierr = IGA_Grid_GetScatterG2L(grid,&iga->g2l);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject)iga->g2l);CHKERRQ(ierr);
     ierr = IGA_Grid_GetScatterL2G(grid,&iga->l2g);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject)iga->l2g);CHKERRQ(ierr);
-    /* build global to natural vector scatter */
-    ierr = IGA_Grid_GetScatterG2N(grid,&iga->g2n);CHKERRQ(ierr);
-    ierr = PetscObjectReference((PetscObject)iga->g2n);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecNatural(grid,VECSTANDARD,&iga->natural);CHKERRQ(ierr);
-    ierr = PetscObjectReference((PetscObject)iga->natural);CHKERRQ(ierr);
+    /* build global <-> natural vector scatter */
+    ierr = IGA_Grid_NewScatterApp(grid,
+                                  iga->geom_sizes,iga->geom_lstart,iga->geom_lwidth,
+                                  &iga->natural,&iga->n2g,&iga->g2n);CHKERRQ(ierr);
     /* destroy the grid context */
     ierr = IGA_Grid_Destroy(&grid);CHKERRQ(ierr);
   }

File src/petigagrid.c

   ierr = AODestroy(&g->aob);CHKERRQ(ierr);
   ierr = ISLocalToGlobalMappingDestroy(&g->lgmap);CHKERRQ(ierr);
   ierr = ISLocalToGlobalMappingDestroy(&g->lgmapb);CHKERRQ(ierr);
+  ierr = VecDestroy(&g->lvec);CHKERRQ(ierr);
+  ierr = VecDestroy(&g->gvec);CHKERRQ(ierr);
   ierr = VecDestroy(&g->nvec);CHKERRQ(ierr);
-  ierr = VecDestroy(&g->gvec);CHKERRQ(ierr);
-  ierr = VecDestroy(&g->lvec);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&g->g2l);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&g->l2g);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&g->g2n);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGA_Grid_GetVecNatural"
-PetscErrorCode IGA_Grid_GetVecNatural(IGA_Grid g,const VecType vtype,Vec *nvec)
+#define __FUNCT__ "IGA_Grid_GetVecLocal"
+PetscErrorCode IGA_Grid_GetVecLocal(IGA_Grid g,const VecType vtype,Vec *lvec)
 {
-  Vec            gvec;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(g,1);
   if (vtype) PetscValidCharPointer(vtype,2);
-  PetscValidPointer(nvec,3);
-  if (!g->nvec) {
-    ierr = IGA_Grid_GetVecGlobal(g,vtype,&gvec);CHKERRQ(ierr);
-    ierr = VecDuplicate(gvec,&g->nvec);CHKERRQ(ierr);
+  PetscValidPointer(lvec,3);
+  if (!g->lvec) {
+    const PetscInt *width = g->ghost_width;
+    PetscInt n  = width[0]*width[1]*width[2];
+    PetscInt bs = g->dof;
+    ierr = VecCreate(PETSC_COMM_SELF,&g->lvec);CHKERRQ(ierr);
+    ierr = VecSetSizes(g->lvec,n*bs,n*bs);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(g->lvec,bs);CHKERRQ(ierr);
+    ierr = VecSetType(g->lvec,vtype?vtype:VECSTANDARD);CHKERRQ(ierr);
   }
-  *nvec = g->nvec;
+  *lvec = g->lvec;
   PetscFunctionReturn(0);
 }
 
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGA_Grid_GetVecLocal"
-PetscErrorCode IGA_Grid_GetVecLocal(IGA_Grid g,const VecType vtype,Vec *lvec)
+#define __FUNCT__ "IGA_Grid_GetVecNatural"
+PetscErrorCode IGA_Grid_GetVecNatural(IGA_Grid g,const VecType vtype,Vec *nvec)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(g,1);
   if (vtype) PetscValidCharPointer(vtype,2);
-  PetscValidPointer(lvec,3);
-  if (!g->lvec) {
-    const PetscInt *width = g->ghost_width;
-    PetscInt n  = width[0]*width[1]*width[2];
-    PetscInt bs = g->dof;
-    ierr = VecCreate(PETSC_COMM_SELF,&g->lvec);CHKERRQ(ierr);
-    ierr = VecSetSizes(g->lvec,n*bs,n*bs);CHKERRQ(ierr);
-    ierr = VecSetBlockSize(g->lvec,bs);CHKERRQ(ierr);
-    ierr = VecSetType(g->lvec,vtype?vtype:VECSTANDARD);CHKERRQ(ierr);
+  PetscValidPointer(nvec,3);
+  if (!g->nvec) {
+    Vec gvec;
+    ierr = IGA_Grid_GetVecGlobal(g,vtype,&gvec);CHKERRQ(ierr);
+    ierr = VecDuplicate(gvec,&g->nvec);CHKERRQ(ierr);
   }
-  *lvec = g->lvec;
+  *nvec = g->nvec;
   PetscFunctionReturn(0);
 }
 
   PetscValidPointer(g,1);
   PetscValidPointer(l2g,2);
   if (!g->l2g) {
+    IS isglobal,islocal;
+    Vec gvec,lvec;
+    PetscInt nlocal,*ilocal,start;
+    /* local non-ghosted grid */
     const PetscInt *lstart = g->local_start;
     const PetscInt *lwidth = g->local_width;
-    const PetscInt *gstart = g->ghost_start;
-    const PetscInt *gwidth = g->ghost_width;
-    /* local non-ghosted grid */
     PetscInt ilstart = lstart[0], ilend = lstart[0]+lwidth[0];
     PetscInt jlstart = lstart[1], jlend = lstart[1]+lwidth[1];
     PetscInt klstart = lstart[2], klend = lstart[2]+lwidth[2];
     /* local ghosted grid */
+    const PetscInt *gstart = g->ghost_start;
+    const PetscInt *gwidth = g->ghost_width;
     PetscInt igstart = gstart[0], igend = gstart[0]+gwidth[0];
     PetscInt jgstart = gstart[1], jgend = gstart[1]+gwidth[1];
     PetscInt kgstart = gstart[2], kgend = gstart[2]+gwidth[2];
-    IS isglobal,islocal;
-    Vec gvec,lvec;
-    PetscInt start,nlocal,*ilocal;
-    PetscInt c,i,j,k,pos = 0,index = 0, bs = g->dof;
+    /* */
+    PetscInt c,i,j,k,pos = 0,index = 0,bs = g->dof;
     ierr = IGA_Grid_GetVecGlobal(g,VECSTANDARD,&gvec);CHKERRQ(ierr);
     ierr = IGA_Grid_GetVecLocal(g,VECSTANDARD,&lvec);CHKERRQ(ierr);
     ierr = VecGetLocalSize(gvec,&nlocal);CHKERRQ(ierr);
     for (k=kgstart; k<kgend; k++)
       for (j=jgstart; j<jgend; j++)
         for (i=igstart; i<igend; i++, index++)
-          if (i>=ilstart && i<ilend && j>=jlstart && j<jlend && k>=klstart && k<klend)
+          if (i>=ilstart && i<ilend &&
+              j>=jlstart && j<jlend &&
+              k>=klstart && k<klend)
             for (c=0; c<bs; c++)
               ilocal[pos++] = c + bs*index;
     ierr = ISCreateGeneral(PETSC_COMM_SELF,nlocal,ilocal,PETSC_OWN_POINTER,&islocal);CHKERRQ(ierr);
   PetscValidPointer(g2n,2);
   if (!g->g2n) {
     IS isnatural,isglobal;
-    Vec gvec;
+    Vec gvec,nvec;
     PetscInt start,nlocal,*inatural;
-    ierr = IGA_Grid_GetVecGlobal(g,VECSTANDARD,&gvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecGlobal (g,VECSTANDARD,&gvec);CHKERRQ(ierr);
+    ierr = IGA_Grid_GetVecNatural(g,VECSTANDARD,&nvec);CHKERRQ(ierr);
     ierr = IGA_Grid_LocalIndices(g,g->dof,&nlocal,&inatural);CHKERRQ(ierr);
     ierr = VecGetOwnershipRange(gvec,&start,PETSC_NULL);CHKERRQ(ierr);
     ierr = ISCreateStride(g->comm,nlocal,start,1,&isglobal);CHKERRQ(ierr);
     ierr = ISCreateGeneral(g->comm,nlocal,inatural,PETSC_OWN_POINTER,&isnatural);CHKERRQ(ierr);
-    ierr = VecScatterCreate(gvec,isglobal,gvec,isnatural,&g->g2n);CHKERRQ(ierr);
+    ierr = VecScatterCreate(gvec,isglobal,nvec,isnatural,&g->g2n);CHKERRQ(ierr);
     ierr = ISDestroy(&isglobal);CHKERRQ(ierr);
     ierr = ISDestroy(&isnatural);CHKERRQ(ierr);
   }
   *g2n = g->g2n;
   PetscFunctionReturn(0);
 }
+
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGA_Grid_NewScatterApp"
+PetscErrorCode IGA_Grid_NewScatterApp(IGA_Grid g,
+                                      const PetscInt sizes[],
+                                      const PetscInt start[],
+                                      const PetscInt width[],
+                                      Vec        *avec,
+                                      VecScatter *a2g,
+                                      VecScatter *g2a)
+{
+  const char*    vtype = VECSTANDARD;
+  Vec            gvec;
+  Vec            nvec;
+  VecScatter     n2g;
+  VecScatter     g2n;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(g,1);
+  PetscValidPointer(avec,2);
+  PetscValidPointer(a2g,3);
+  PetscValidPointer(g2a,4);
+  /* global vector */
+  {
+    ierr = IGA_Grid_GetVecGlobal(g,vtype,&gvec);CHKERRQ(ierr);
+    ierr = VecGetType(gvec,&vtype);CHKERRQ(ierr);
+  }
+  /* natural vector */
+  {
+    PetscInt n  = width[0]*width[1]*width[2];
+    PetscInt N  = sizes[0]*sizes[1]*sizes[2];
+    PetscInt bs = g->dof;
+    ierr = VecCreate(g->comm,&nvec);CHKERRQ(ierr);
+    ierr = VecSetSizes(nvec,n*bs,N*bs);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(nvec,bs);CHKERRQ(ierr);
+    ierr = VecSetType(nvec,vtype);CHKERRQ(ierr);
+  }
+  /* natural -> global scatter */
+  {
+    IS isnatural,isglobal;
+    PetscInt nlocal,*inatural,nstart;
+    /* global grid strides */
+    PetscInt jstride = sizes[0];
+    PetscInt kstride = sizes[0]*sizes[1];
+    /* local non-ghosted grid */
+    const PetscInt *lstart = g->local_start;
+    const PetscInt *lwidth = g->local_width;
+    PetscInt ilstart = lstart[0], ilend = lstart[0]+lwidth[0];
+    PetscInt jlstart = lstart[1], jlend = lstart[1]+lwidth[1];
+    PetscInt klstart = lstart[2], klend = lstart[2]+lwidth[2];
+    /* */
+    PetscInt c,i,j,k,pos = 0,bs = g->dof;
+    ierr = VecGetLocalSize(gvec,&nlocal);CHKERRQ(ierr);
+    ierr = VecGetOwnershipRange(gvec,&nstart,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscMalloc(nlocal*sizeof(PetscInt),&inatural);CHKERRQ(ierr);
+    for (k=klstart; k<klend; k++)
+      for (j=jlstart; j<jlend; j++)
+        for (i=ilstart; i<ilend; i++)
+          for (c=0; c<bs; c++)
+            inatural[pos++] = c + bs*(i + j * jstride + k * kstride);
+    ierr = ISCreateGeneral(g->comm,nlocal,inatural,PETSC_OWN_POINTER,&isnatural);CHKERRQ(ierr);
+    ierr = ISCreateStride(g->comm,nlocal,nstart,1,&isglobal);CHKERRQ(ierr);
+    ierr = VecScatterCreate(nvec,isnatural,gvec,isglobal,&n2g);CHKERRQ(ierr);
+    ierr = ISDestroy(&isnatural);CHKERRQ(ierr);
+    ierr = ISDestroy(&isglobal);CHKERRQ(ierr);
+  }
+  /* global -> natural scatter */
+  {
+    IS isglobal,isnatural;
+    PetscInt nlocal,*iglobal,nstart;
+    /* local non-ghosted grid */
+    const PetscInt *lstart = start;
+    const PetscInt *lwidth = width;
+    PetscInt ilstart = lstart[0], ilend = lstart[0]+lwidth[0];
+    PetscInt jlstart = lstart[1], jlend = lstart[1]+lwidth[1];
+    PetscInt klstart = lstart[2], klend = lstart[2]+lwidth[2];
+    /* local ghosted grid */
+    const PetscInt *gstart = g->ghost_start;
+    const PetscInt *gwidth = g->ghost_width;
+    PetscInt igstart = gstart[0], igend = gstart[0]+gwidth[0];
+    PetscInt jgstart = gstart[1], jgend = gstart[1]+gwidth[1];
+    PetscInt kgstart = gstart[2], kgend = gstart[2]+gwidth[2];
+    /* */
+    PetscInt c,i,j,k,pos = 0,index = 0,bs = g->dof;
+    ierr = VecGetLocalSize(nvec,&nlocal);CHKERRQ(ierr);
+    ierr = VecGetOwnershipRange(nvec,&nstart,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscMalloc(nlocal*sizeof(PetscInt),&iglobal);CHKERRQ(ierr);
+    for (k=kgstart; k<kgend; k++)
+      for (j=jgstart; j<jgend; j++)
+        for (i=igstart; i<igend; i++, index++)
+          if (i>=ilstart && i<ilend &&
+              j>=jlstart && j<jlend &&
+              k>=klstart && k<klend)
+            for (c=0; c<bs; c++)
+              iglobal[pos++] = c + bs*index;
+    ierr = IGA_Grid_GetLGMap(g,&g->lgmap);CHKERRQ(ierr);
+    ierr = ISLocalToGlobalMappingApply(g->lgmap,nlocal,iglobal,iglobal);CHKERRQ(ierr);
+    ierr = ISCreateGeneral(g->comm,nlocal,iglobal,PETSC_OWN_POINTER,&isglobal);CHKERRQ(ierr);
+    ierr = ISCreateStride(g->comm,nlocal,nstart,1,&isnatural);CHKERRQ(ierr);
+    ierr = VecScatterCreate(gvec,isglobal,nvec,isnatural,&g2n);CHKERRQ(ierr);
+    ierr = ISDestroy(&isnatural);CHKERRQ(ierr);
+    ierr = ISDestroy(&isglobal);CHKERRQ(ierr);
+  }
+  *avec = nvec;
+  *a2g  = n2g;
+  *g2a  = g2n;
+  PetscFunctionReturn(0);
+}

File src/petigagrid.h

   PetscInt   ghost_width[3];
   AO         ao,aob;
   LGMap      lgmap,lgmapb;
-  Vec        nvec,gvec,lvec;
+  Vec        lvec,gvec,nvec;
   VecScatter g2l,l2g,g2n;
 };
 
 PETSC_EXTERN PetscErrorCode IGA_Grid_SetLGMapBlock(IGA_Grid,LGMap);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetLGMapBlock(IGA_Grid,LGMap*);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetLGMap(IGA_Grid,LGMap*);
+PETSC_EXTERN PetscErrorCode IGA_Grid_GetVecLocal  (IGA_Grid,const VecType,Vec*);
+PETSC_EXTERN PetscErrorCode IGA_Grid_GetVecGlobal (IGA_Grid,const VecType,Vec*);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetVecNatural(IGA_Grid,const VecType,Vec*);
-PETSC_EXTERN PetscErrorCode IGA_Grid_GetVecGlobal (IGA_Grid,const VecType,Vec*);
-PETSC_EXTERN PetscErrorCode IGA_Grid_GetVecLocal  (IGA_Grid,const VecType,Vec*);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetScatterG2L(IGA_Grid,VecScatter*);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetScatterL2G(IGA_Grid,VecScatter*);
 PETSC_EXTERN PetscErrorCode IGA_Grid_GetScatterG2N(IGA_Grid,VecScatter*);
 
+PETSC_EXTERN PetscErrorCode IGA_Grid_NewScatterApp(IGA_Grid g,
+                                                   const PetscInt[],
+                                                   const PetscInt[],
+                                                   const PetscInt[],
+                                                   Vec*,VecScatter*,VecScatter*);
+
 #if PETSC_VERSION_(3,2,0)
 PETSC_EXTERN_CXX_BEGIN
 #endif

File src/petigaio.c

 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGA_NewGridIO"
-static PetscErrorCode IGA_NewGridIO(IGA iga,PetscInt bs,IGA_Grid *grid)
+#define __FUNCT__ "IGA_NewGridGeom"
+static PetscErrorCode IGA_NewGridGeom(IGA iga,IGA_Grid *grid)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(grid,3);
   {
     MPI_Comm comm    = ((PetscObject)iga)->comm;
+    PetscInt dim     = iga->dim;
     PetscInt *sizes  = iga->geom_sizes;
     PetscInt *lstart = iga->geom_lstart;
     PetscInt *lwidth = iga->geom_lwidth;
     PetscInt *gstart = iga->geom_gstart;
     PetscInt *gwidth = iga->geom_gwidth;
     ierr = IGA_Grid_Create(comm,grid);CHKERRQ(ierr);
-    ierr = IGA_Grid_Init(*grid,iga->dim,bs,sizes,lstart,lwidth,gstart,gwidth);CHKERRQ(ierr);
+    ierr = IGA_Grid_Init(*grid,dim,dim+1,sizes,lstart,lwidth,gstart,gwidth);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
   if (dim < 1)
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,
             "Must call IGASetSpatialDim() first");
+
+  iga->geometry = PETSC_FALSE;
+  iga->rational = PETSC_FALSE;
+  ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
+  ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
+  ierr = VecDestroy(&iga->geom_vec);CHKERRQ(ierr);
+
   {
     IGA_Grid grid;
-    ierr = IGA_NewGridIO(iga,dim+1,&grid);CHKERRQ(ierr);
+    ierr = IGA_NewGridGeom(iga,&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_Destroy(&grid);CHKERRQ(ierr);
   }
   ierr = PetscObjectReference((PetscObject)lvec);CHKERRQ(ierr);
-  ierr = VecDestroy(&iga->geom_vec);CHKERRQ(ierr);
   iga->geom_vec = lvec;
 
   /* viewer -> natural*/
     ierr = VecGetSize(iga->geom_vec,&n);CHKERRQ(ierr);
     ierr = VecGetBlockSize(iga->geom_vec,&bs);CHKERRQ(ierr);
     nnp = n / bs; dim = bs - 1;
-    ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
-    ierr = PetscFree(iga->geometryW);CHKERRQ(ierr);
     ierr = PetscMalloc1(nnp*dim,PetscReal,&iga->geometryX);CHKERRQ(ierr);
     ierr = PetscMalloc1(nnp,    PetscReal,&iga->geometryW);CHKERRQ(ierr);
     X = iga->geometryX; W = iga->geometryW;
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,
             "Must call IGASetSpatialDim() first");
   {
-    IGA_Grid   grid;
-    ierr = IGA_NewGridIO(iga,dim+1,&grid);CHKERRQ(ierr);
+    IGA_Grid grid;
+    ierr = IGA_NewGridGeom(iga,&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);
   /* 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 -> naural */
+  /* 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 */
   PetscFunctionReturn(0);
 }
 
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGA_AnyPeriodic"
-static PetscBool IGA_AnyPeriodic(IGA iga)
-{
-  PetscInt i, dim = iga->dim;
-  for (i=0; i<dim; i++)
-    if (iga->axis[i]->periodic)
-      return PETSC_TRUE;
-  return PETSC_FALSE;
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGA_GridIO_destroy"
-static PetscErrorCode IGA_GridIO_destroy(void *ptr)
-{
-  IGA_Grid       grid = (IGA_Grid)ptr;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(grid,1);
-  ierr = IGA_Grid_Destroy(&grid);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGA_GetGridIO_Vec"
-static PetscErrorCode IGA_GetGridIO_Vec(IGA iga,IGA_Grid *gridio)
-{
-  const char     key[] = "_GridIO_Vec";
-  MPI_Comm       comm;
-  IGA_Grid       grid;
-  PetscContainer container,tmp;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidPointer(gridio,4);
-  ierr = PetscObjectQuery((PetscObject)iga,key,(PetscObject*)&container);CHKERRQ(ierr);
-  if (!container) {
-    ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
-    ierr = PetscContainerCreate(comm,&container);CHKERRQ(ierr);
-    ierr = PetscObjectCompose((PetscObject)iga,key,(PetscObject)container);CHKERRQ(ierr);
-    tmp = container;
-    ierr = PetscContainerDestroy(&tmp);CHKERRQ(ierr);
-  }
-  ierr = PetscContainerGetPointer(container,(void**)&grid);CHKERRQ(ierr);
-  if (!grid) {
-    ierr = IGA_NewGridIO(iga,iga->dof,&grid);CHKERRQ(ierr);
-    ierr = PetscContainerSetPointer(container,(void*)grid);CHKERRQ(ierr);
-    ierr = PetscContainerSetUserDestroy(container,IGA_GridIO_destroy);;CHKERRQ(ierr);
-  }
-  *gridio = grid;
-  PetscFunctionReturn(0);
-}
-
 #undef  __FUNCT__
 #define __FUNCT__ "IGALoadVec"
 PetscErrorCode IGALoadVec(IGA iga,Vec vec,PetscViewer viewer)
 {
-  Vec            nvec;
+  Vec            natural;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscCheckSameComm(iga,1,viewer,3);
   IGACheckSetUp(iga,1);
 
-  if (!IGA_AnyPeriodic(iga)) {
-    ierr = IGAGetNaturalVec(iga,&nvec);
-    ierr = VecLoad(nvec,viewer);CHKERRQ(ierr);
-    ierr = IGANaturalToGlobal(iga,nvec,vec);
-  } else {
-    IGA_Grid   grid;
-    Vec        gvec,lvec;
-    VecScatter g2n,g2l;
-    ierr = IGA_GetGridIO_Vec(iga,&grid);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecNatural(grid,VECSTANDARD,&nvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecGlobal (grid,VECSTANDARD,&gvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecLocal  (grid,VECSTANDARD,&lvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetScatterG2N(grid,&g2n);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetScatterG2L(grid,&g2l);CHKERRQ(ierr);
-    ierr = VecLoad(nvec,viewer);CHKERRQ(ierr);
-    ierr = VecScatterBegin(g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterBegin(g2l,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (g2l,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = IGALocalToGlobal(iga,lvec,vec,INSERT_VALUES);CHKERRQ(ierr);
-  }
+  ierr = IGAGetNaturalVec(iga,&natural);
+  ierr = VecLoad(natural,viewer);CHKERRQ(ierr);
+  ierr = IGANaturalToGlobal(iga,natural,vec);
+
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "IGASaveVec"
 PetscErrorCode IGASaveVec(IGA iga,Vec vec,PetscViewer viewer)
 {
-  Vec            nvec;
+  Vec            natural;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscCheckSameComm(iga,1,viewer,3);
   IGACheckSetUp(iga,1);
 
-  if (!IGA_AnyPeriodic(iga)) {
-    ierr = IGAGetNaturalVec(iga,&nvec);
-    ierr = IGAGlobalToNatural(iga,vec,nvec);
-  } else {
-    IGA_Grid   grid;
-    Vec        gvec,lvec;
-    VecScatter g2n,l2g;
-    ierr = IGA_GetGridIO_Vec(iga,&grid);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecLocal  (grid,VECSTANDARD,&lvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecGlobal (grid,VECSTANDARD,&gvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetVecNatural(grid,VECSTANDARD,&nvec);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetScatterL2G(grid,&l2g);CHKERRQ(ierr);
-    ierr = IGA_Grid_GetScatterG2N(grid,&g2n);CHKERRQ(ierr);
-    ierr = IGAGlobalToLocal(iga,vec,lvec);CHKERRQ(ierr);
-    ierr = VecScatterBegin(l2g,lvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (l2g,lvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterBegin(g2n,gvec,nvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (g2n,gvec,nvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  }
-  ierr = PetscObjectSetName((PetscObject)nvec,((PetscObject)vec)->name);CHKERRQ(ierr);
-  ierr = VecView(nvec,viewer);CHKERRQ(ierr);
+  ierr = IGAGetNaturalVec(iga,&natural);
+  ierr = IGAGlobalToNatural(iga,vec,natural);
+  ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)vec)->name);CHKERRQ(ierr);
+  ierr = VecView(natural,viewer);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }

File src/petigavec.c

   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(nvec,VEC_CLASSID,2);
   PetscValidHeaderSpecific(gvec,VEC_CLASSID,3);
-  ierr = VecScatterBegin(iga->g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (iga->g2n,nvec,gvec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  ierr = VecScatterBegin(iga->n2g,nvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd  (iga->n2g,nvec,gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }