Always use a taggable set (not root) as the domain.

I'm not sure if this is just a MOAB quirk, but you can't set a tag on the root set.
jedbrown · Nov 29, 2009 · 31f4feb · 31f4feb
1 parent ca14d80
commit 31f4feb
Show file tree

Hide file tree

Showing 8 changed files with 44 additions and 11 deletions.
diff --git a/include/dohpmesh.h b/include/dohpmesh.h
@@ -105,6 +105,7 @@ EXTERN dErr dMeshGetTagName(dMesh m,dMeshTag tag,char **name);
 EXTERN dErr dMeshLoad(dMesh m);
 EXTERN dErr dMeshSetInFile(dMesh,const char fname[],const char opt[]);
 EXTERN dErr dMeshGetRoot(dMesh mesh,dMeshESH *inroot);
+EXTERN dErr dMeshSetDuplicateEntsOnly(dMesh mesh,dMeshESH set,dMeshESH *copy);
 EXTERN dErr dMeshCreate(MPI_Comm comm,dMesh *inm);
 EXTERN dErr dMeshDestroy(dMesh);
 EXTERN dErr dMeshView(dMesh,PetscViewer);

diff --git a/src/fs/mesh/interface/mesh.c b/src/fs/mesh/interface/mesh.c
@@ -313,6 +313,28 @@ dErr dMeshGetRoot(dMesh mesh,dMeshESH *inroot)
   dFunctionReturn(0);
 }
 
+/* Makes a new set containing all entities in original, but not subsets */
+dErr dMeshSetDuplicateEntsOnly(dMesh mesh,dMeshESH set,dMeshESH *copy)
+{
+  iMesh_Instance mi = mesh->mi;
+  dIInt   ierr,islist;
+  dInt    size;
+  dMeshEH *ents;
+  dErr    err;
+
+  dFunctionBegin;
+  dValidHeader(mesh,dMESH_COOKIE,1);
+  dValidPointer(copy,3);
+  err = dMeshGetNumEnts(mesh,set,dTYPE_ALL,dTOPO_ALL,&size);dCHK(err);
+  err = dMallocA(size,&ents);dCHK(err);
+  err = dMeshGetEnts(mesh,set,dTYPE_ALL,dTOPO_ALL,ents,size,NULL);dCHK(err);
+  iMesh_isList(mi,set,&islist,&ierr);dICHK(mi,ierr);
+  err = dMeshSetCreate(mesh,islist?dMESHSET_ORDERED:dMESHSET_UNORDERED,copy);dCHK(err);
+  err = dMeshSetAddEnts(mesh,*copy,ents,size);dCHK(err);
+  err = dFree(ents);dCHK(err);
+  dFunctionReturn(0);
+}
+
 /**
 * This function allocates memory for the tag name.  It should be freed with PetscFree()
 */

diff --git a/src/fs/tests/bu.c b/src/fs/tests/bu.c
@@ -58,6 +58,7 @@ static dErr BUSetFromOptions(BU bu)
   err = dMeshLoad(mesh);dCHK(err);dCHK(err);
   bu->mesh = mesh;
   err = dMeshGetRoot(mesh,&domain);dCHK(err);
+  err = dMeshSetDuplicateEntsOnly(mesh,domain,&domain);dCHK(err);
 
   err = dJacobiCreate(bu->comm,&jac);dCHK(err);
   err = dJacobiSetDegrees(jac,9,2);dCHK(err);
@@ -69,7 +70,7 @@ static dErr BUSetFromOptions(BU bu)
   err = dMeshCreateRuleTagIsotropic(mesh,domain,jac,"bu_efs_degree",bu->constBDeg,&dtag);dCHK(err);
 
   err = dFSCreate(bu->comm,&fs);dCHK(err);
-  err = dFSSetMesh(fs,mesh,0);dCHK(err);
+  err = dFSSetMesh(fs,mesh,domain);dCHK(err);
   err = dFSSetRuleTag(fs,jac,rtag);dCHK(err);
   err = dFSSetDegree(fs,jac,dtag);dCHK(err);
   bu->fs = fs;

diff --git a/src/fs/tests/elast.c b/src/fs/tests/elast.c
@@ -167,7 +167,8 @@ static dErr ElastSetFromOptions(Elast elt)
   err = dMeshSetFromOptions(mesh);dCHK(err);
   err = dMeshLoad(mesh);dCHK(err);dCHK(err);
   elt->mesh = mesh;
-  domain = 0;                   /* Root set in MOAB */
+  err = dMeshGetRoot(mesh,&domain);dCHK(err); /* Need a taggable set */
+  err = dMeshSetDuplicateEntsOnly(mesh,domain,&domain);dCHK(err);
 
   err = dJacobiCreate(elt->comm,&jac);dCHK(err);
   err = dJacobiSetDegrees(jac,9,2);dCHK(err);
@@ -180,7 +181,7 @@ static dErr ElastSetFromOptions(Elast elt)
 
   err = dFSCreate(elt->comm,&fs);dCHK(err);
   err = dFSSetBlockSize(fs,3);dCHK(err);
-  err = dFSSetMesh(fs,mesh,0);dCHK(err);
+  err = dFSSetMesh(fs,mesh,domain);dCHK(err);
   err = dFSSetRuleTag(fs,jac,rtag);dCHK(err);
   err = dFSSetDegree(fs,jac,dtag);dCHK(err);
   err = dFSRegisterBoundary(fs,100,dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);

diff --git a/src/fs/tests/ellip.c b/src/fs/tests/ellip.c
@@ -221,7 +221,8 @@ static dErr EllipSetFromOptions(Ellip elp)
     err = dMeshMorph(mesh,&Morph,ctx);dCHK(err);
   }
   elp->mesh = mesh;
-  domain = 0;                   /* Root set in MOAB */
+  err = dMeshGetRoot(mesh,&domain);dCHK(err); /* Need a taggable set */
+  err = dMeshSetDuplicateEntsOnly(mesh,domain,&domain);dCHK(err);
 
   err = dJacobiCreate(elp->comm,&jac);dCHK(err);
   err = dJacobiSetDegrees(jac,9,2);dCHK(err);
@@ -233,7 +234,7 @@ static dErr EllipSetFromOptions(Ellip elp)
   err = dMeshCreateRuleTagIsotropic(mesh,domain,jac,"ellip_efs_degree",elp->constBDeg,&dtag);dCHK(err);
 
   err = dFSCreate(elp->comm,&fs);dCHK(err);
-  err = dFSSetMesh(fs,mesh,0);dCHK(err);
+  err = dFSSetMesh(fs,mesh,domain);dCHK(err);
   err = dFSSetRuleTag(fs,jac,rtag);dCHK(err);
   err = dFSSetDegree(fs,jac,dtag);dCHK(err);
   err = dFSRegisterBoundary(fs,100,dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);

diff --git a/src/fs/tests/ex1.c b/src/fs/tests/ex1.c
@@ -137,7 +137,7 @@ static dErr tagHexes(dMesh mesh,dMeshTag *intag)
   if (gopt.constBDeg) {
     for (dInt i=0; i<ALEN(hexDegree); i++) hexDegree[i] = gopt.constBDeg;
   }
-  err = dMeshTagCreateTemp(mesh,"anisotropic",3,dDATA_INT,&tag);dCHK(err);
+  err = dMeshTagCreate(mesh,"ex1_anisotropic",3,dDATA_INT,&tag);dCHK(err);
   /* tag edges and faces with high values, currently needed to propogate degrees (will overwrite high values) */
   for (dEntType type=dTYPE_VERTEX; type<=dTYPE_FACE; type++) {
     err = dMeshGetEnts(mesh,0,type,dTOPO_ALL,ents,ALEN(ents),&nents);
@@ -509,7 +509,10 @@ int main(int argc,char *argv[])
   err = createHexMesh(mi);dCHK(err);
   if (showmesh) {err = dMeshView(mesh,viewer);dCHK(err);}
   err = verifyAdjacencies(mesh);dCHK(err);
-  iMesh_getRootSet(mi,&domain,&err);dICHK(mi,err);
+
+  /* The root set cannot be tagged (at least in MOAB) */
+  err = dMeshGetRoot(mesh,&domain);dCHK(err);
+  err = dMeshSetDuplicateEntsOnly(mesh,domain,&domain);dCHK(err);
 
   err = dJacobiCreate(comm,&jac);dCHK(err);
   err = dJacobiSetDegrees(jac,15,4);dCHK(err);

diff --git a/src/fs/tests/fsex4.c b/src/fs/tests/fsex4.c
@@ -95,6 +95,7 @@ int main(int argc,char *argv[])
   dViewer  viewer;
   dJacobi  jac;
   dMeshTag rtag,dtag;
+  dMeshESH active;
   dTruth   read;
   Vec      X;
 
@@ -115,11 +116,13 @@ int main(int argc,char *argv[])
   err = dJacobiSetDegrees(jac,9,2);dCHK(err);
   err = dJacobiSetFromOptions(jac);dCHK(err);
   err = dJacobiSetUp(jac);dCHK(err);
+  err = dMeshGetRoot(mesh,&active);dCHK(err); /* Need a taggable set */
+  err = dMeshSetDuplicateEntsOnly(mesh,active,&active);dCHK(err);
 
-  err = dMeshCreateRuleTagIsotropic(mesh,0,jac,"fsex4_rule_degree",ex4->rdeg,&rtag);dCHK(err);
-  err = dMeshCreateRuleTagIsotropic(mesh,0,jac,"fsex4_efs_degree",ex4->bdeg,&dtag);dCHK(err);
+  err = dMeshCreateRuleTagIsotropic(mesh,active,jac,"fsex4_rule_degree",ex4->rdeg,&rtag);dCHK(err);
+  err = dMeshCreateRuleTagIsotropic(mesh,active,jac,"fsex4_efs_degree",ex4->bdeg,&dtag);dCHK(err);
   err = dFSCreate(comm,&fs);dCHK(err);
-  err = dFSSetMesh(fs,mesh,0);dCHK(err);
+  err = dFSSetMesh(fs,mesh,active);dCHK(err);
   err = dFSSetRuleTag(fs,jac,rtag);dCHK(err);
   err = dFSSetDegree(fs,jac,dtag);dCHK(err);
   err = dFSSetFromOptions(fs);dCHK(err);

diff --git a/src/fs/tests/stokes.c b/src/fs/tests/stokes.c
@@ -257,7 +257,8 @@ static dErr StokesSetFromOptions(Stokes stk)
   err = dMeshSetFromOptions(mesh);dCHK(err);
   err = dMeshLoad(mesh);dCHK(err);dCHK(err);
   stk->mesh = mesh;
-  err = dMeshGetRoot(mesh,&domain);dCHK(err);
+  err = dMeshGetRoot(mesh,&domain);dCHK(err); /* Need a taggable set */
+  err = dMeshSetDuplicateEntsOnly(mesh,domain,&domain);dCHK(err);
 
   err = dJacobiCreate(stk->comm,&jac);dCHK(err);
   err = dJacobiSetDegrees(jac,9,2);dCHK(err);