Commits

Vijay Mahadevan  committed 9697b56

Adding several important public API functions to DMMoab.

1) DMMoabSetGlobalFieldVector - Use when the number of fields > 1 and want to assign a single global solution vector to be split in terms of tag arrays.
2) DMMoabGetDofs - Local/Blocked variants to get the individual field-independent vertex GID numbering. This is useful for traversals based on nodes directly. Used in FD-type computations.
3) Some minor bug fixes.

  • Participants
  • Parent commits a3d2f62

Comments (0)

Files changed (2)

File include/petscdmmoab.h

 PETSC_EXTERN PetscErrorCode DMMoabCreateMatrix(DM dm, MatType mtype,Mat *J);
 PETSC_EXTERN PetscErrorCode DMMoabGetVecTag(Vec vec,moab::Tag *tag);
 PETSC_EXTERN PetscErrorCode DMMoabGetVecRange(Vec vec,moab::Range *range);
+
 PETSC_EXTERN PetscErrorCode DMMoabSetFieldVector(DM, PetscInt, Vec);
+PETSC_EXTERN PetscErrorCode DMMoabSetGlobalFieldVector(DM, Vec);
 
 PETSC_EXTERN PetscErrorCode DMMoabSetFields(DM dm,PetscInt nfields,const char** fields);
 PETSC_EXTERN PetscErrorCode DMMoabGetFieldDof(DM dm,moab::EntityHandle point, PetscInt field,PetscInt* dof);
 PETSC_EXTERN PetscErrorCode DMMoabGetFieldDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points, PetscInt field,PetscInt* dof);
 PETSC_EXTERN PetscErrorCode DMMoabGetLocalFieldDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points, PetscInt field,PetscInt* dof);
+PETSC_EXTERN PetscErrorCode DMMoabGetDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof);
+PETSC_EXTERN PetscErrorCode DMMoabGetLocalDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof);
+PETSC_EXTERN PetscErrorCode DMMoabGetDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof);
+PETSC_EXTERN PetscErrorCode DMMoabGetLocalDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof);
 
 /* discretization and assembly specific DMMoab interface functions */
 PETSC_EXTERN PetscErrorCode DMMoabGetElementConnectivity(DM dm,moab::EntityHandle ehandle,PetscInt* nconn,const moab::EntityHandle **conn);

File src/dm/impls/moab/dmmoab.cxx

 {
   DM_Moab        *dmmoab;
   moab::Tag     vtag,ntag;
-  PetscScalar   *varray;
+  const PetscScalar *varray;
+  PetscScalar *farray;
   moab::ErrorCode merr;
   PetscErrorCode  ierr;
-  PetscSection section;
   PetscInt doff;
   std::string tag_name;
   moab::Range::iterator iter;
   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
   dmmoab = (DM_Moab*)(dm)->data;
 
-  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-
   /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
   merr = dmmoab->mbiface->tag_get_handle(dmmoab->fields[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
                                           moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT);MBERRNM(merr);
 
   merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
   if (!tag_name.length() && merr !=moab::MB_SUCCESS) {
-    ierr = DMMoabVecGetArray(dm,fvec,&varray);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(fvec,&varray);CHKERRQ(ierr);
+
     for(iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++) {
       moab::EntityHandle vtx = (*iter);
 
       /* get field dof index */
-      ierr = PetscSectionGetFieldOffset(section, vtx, ifield, &doff);
+      ierr = DMMoabGetFieldDof(dm, vtx, ifield, &doff);
 
       /* use the entity handle and the Dof index to set the right value */
       merr = dmmoab->mbiface->tag_set_data(ntag, &vtx, 1, (const void*)&varray[doff]);MBERRNM(merr);
     }
-    ierr = DMMoabVecRestoreArray(dm,fvec,&varray);CHKERRQ(ierr);
+    ierr = VecRestoreArrayRead(fvec,&varray);CHKERRQ(ierr);
   }
   else {
-    ierr = PetscMalloc(dmmoab->nloc*sizeof(PetscScalar),&varray);CHKERRQ(ierr);
+    ierr = PetscMalloc(dmmoab->nloc*sizeof(PetscScalar),&farray);CHKERRQ(ierr);
     /* we are using a MOAB Vec - directly copy the tag data to new one */
-    merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray);MBERRNM(merr);
-    merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)varray);MBERRNM(merr);
+    merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)farray);MBERRNM(merr);
+    merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr);
     /* make sure the parallel exchange for ghosts are done appropriately */
     merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);MBERRNM(merr);
+    ierr = PetscFree(farray);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetGlobalFieldVector"
+PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
+{
+  DM_Moab        *dmmoab;
+  moab::Tag     vtag,ntag;
+  const PetscScalar   *varray;
+  PetscScalar   *farray;
+  moab::ErrorCode merr;
+  PetscErrorCode  ierr;
+  PetscSection section;
+  PetscInt i,doff,ifield;
+  std::string tag_name;
+  moab::Range::iterator iter;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  dmmoab = (DM_Moab*)(dm)->data;
+
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+
+  /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
+  ierr = DMMoabGetVecTag(fvec,&vtag);CHKERRQ(ierr);
+  merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
+  if (!tag_name.length() && merr !=moab::MB_SUCCESS) {
+    /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
+
+    ierr = VecGetArrayRead(fvec,&varray);CHKERRQ(ierr);
+    for (ifield=0; ifield<dmmoab->nfields; ++ifield) {
+
+      /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
+      merr = dmmoab->mbiface->tag_get_handle(dmmoab->fields[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
+                                            moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT);MBERRNM(merr);
+
+      for(iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++) {
+        moab::EntityHandle vtx = (*iter);
+
+        /* get field dof index */
+        ierr = DMMoabGetFieldDof(dm, vtx, ifield, &doff);
+
+        /* use the entity handle and the Dof index to set the right value */
+        merr = dmmoab->mbiface->tag_set_data(ntag, &vtx, 1, (const void*)&varray[doff]);MBERRNM(merr);
+      }
+    }
+    ierr = VecRestoreArrayRead(fvec,&varray);CHKERRQ(ierr);
+  }
+  else {
+    ierr = PetscMalloc(dmmoab->nloc*sizeof(PetscScalar),&farray);CHKERRQ(ierr);
+    ierr = PetscMalloc(dmmoab->nloc*dmmoab->bs*sizeof(PetscScalar),&varray);CHKERRQ(ierr);
+
+    /* we are using a MOAB Vec - directly copy the tag data to new one */
+    merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray);MBERRNM(merr);
+    for (ifield=0; ifield<dmmoab->nfields; ++ifield) {
+
+      /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
+      merr = dmmoab->mbiface->tag_get_handle(dmmoab->fields[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
+                                            moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT);MBERRNM(merr);
+
+      /* we are using a MOAB Vec - directly copy the tag data to new one */
+      for(i=0; i < dmmoab->nloc; i++) {
+        farray[i] = varray[i*dmmoab->bs+ifield];
+      }
+
+      merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr);
+      /* make sure the parallel exchange for ghosts are done appropriately */
+      merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);MBERRNM(merr);
+    }
+    ierr = PetscFree(farray);CHKERRQ(ierr);
     ierr = PetscFree(varray);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 
 
+
 #undef __FUNCT__
 #define __FUNCT__ "DMMoabGetVertexCoordinates"
 PetscErrorCode DMMoabGetVertexCoordinates(DM dm,PetscInt nconn,const moab::EntityHandle *conn,PetscScalar *vpos)
 
 
 #undef __FUNCT__
+#define __FUNCT__ "DMMoabGetDofs"
+PetscErrorCode DMMoabGetDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
+{
+  PetscInt i,f,gid;
+  PetscSection section;
+  PetscErrorCode  ierr;
+  moab::ErrorCode merr;
+  DM_Moab        *dmmoab;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  dmmoab = (DM_Moab*)(dm)->data;
+
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  if (!dof) {
+    ierr = PetscMalloc(sizeof(PetscScalar)*dmmoab->nfields*npoints, &dof);CHKERRQ(ierr);
+  }
+
+  /* first get the local indices */
+  merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,points,npoints,dof);MBERRNM(merr);
+
+  for (i=0; i<npoints; ++i) {
+    gid=dof[i];
+    for (f=0; f<dmmoab->nfields; ++f) {
+      ierr = PetscSectionGetFieldDof(section, gid, f, &dof[i*dmmoab->nfields+f]);CHKERRQ(ierr);
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetLocalDofs"
+PetscErrorCode DMMoabGetLocalDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
+{
+  PetscInt        i,f,offset;
+  PetscErrorCode  ierr;
+  DM_Moab        *dmmoab;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  dmmoab = (DM_Moab*)(dm)->data;
+
+  if (!dof) {
+    ierr = PetscMalloc(sizeof(PetscScalar)*dmmoab->nfields*npoints, &dof);CHKERRQ(ierr);
+  }
+
+  if (dmmoab->bs > 1) {
+    for (f=0; f<dmmoab->nfields; ++f)
+      for (i=0; i<npoints; ++i)
+        dof[i*dmmoab->nfields+f] = (points[i]-1)*dmmoab->bs+f;
+  }
+  else {
+    for (f=0; f<dmmoab->nfields; ++f) {
+      offset = f*dmmoab->n;     /* assume all fields have equal distribution - say all vertex based */
+      for (i=0; i<npoints; ++i)
+        dof[i*dmmoab->nfields+f] = offset+points[i]-1;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetDofsBlocked"
+PetscErrorCode DMMoabGetDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
+{
+  PetscInt i,gid,dofindx;
+  PetscSection section;
+  PetscErrorCode  ierr;
+  moab::ErrorCode merr;
+  DM_Moab        *dmmoab;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  dmmoab = (DM_Moab*)(dm)->data;
+
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  if (!dof) {
+    ierr = PetscMalloc(sizeof(PetscScalar)*npoints, &dof);CHKERRQ(ierr);
+  }
+
+  /* first get the local indices */
+  merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,points,npoints,dof);MBERRNM(merr);
+
+  for (i=0; i<npoints; ++i) {
+    gid=dof[i];
+    ierr = PetscSectionGetFieldDof(section, gid, 0, &dofindx);CHKERRQ(ierr);
+    if (dmmoab->bs > 1)  dof[i]=dofindx/dmmoab->bs;
+    else dof[i]=dofindx;
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetLocalDofsBlocked"
+PetscErrorCode DMMoabGetLocalDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
+{
+  PetscInt        i;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+
+  if (!dof) {
+    ierr = PetscMalloc(sizeof(PetscScalar)*npoints, &dof);CHKERRQ(ierr);
+  }
+
+  for (i=0; i<npoints; ++i)
+    dof[i] = points[i]-1;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMMoab_GetWriteOptions_Private"
 PetscErrorCode DMMoab_GetWriteOptions_Private(PetscInt fsetid, PetscInt numproc, PetscInt dim, MoabWriteMode mode, PetscInt dbglevel, const char* extra_opts, const char** write_opts)
 {