Propogate tensor degree fields down to lower-order entities.

Signed-off-by: Jed Brown <jed@59A2.org>

CMakeLists.txt

@@ -35,7 +35,7 @@ endif (Dohp_USE_DEBUG)
 
 if (PEDANTIC_WARNINGS)
   set (DEFAULT_PEDANTIC_FLAGS "-pedantic -Wall -Wextra -Winline -Wshadow -Wconversion -Wlogical-op -Wmissing-prototypes -Wvla")
-  set (DEFAULT_PEDANTIC_FLAGS "${DEFAULT_PEDANTIC_FLAGS} -Wno-sign-conversion -pedantic -Wstrict-aliasing -fstrict-aliasing")
+  set (DEFAULT_PEDANTIC_FLAGS "${DEFAULT_PEDANTIC_FLAGS} -Wno-sign-conversion -Wwrite-strings -Wstrict-aliasing -Wcast-align -fstrict-aliasing")
   set (DEFAULT_PEDANTIC_FLAGS "${DEFAULT_PEDANTIC_FLAGS} -Wdisabled-optimization -funit-at-a-time")
   #set (DEFAULT_PEDANTIC_FLAGS "${DEFAULT_PEDANTIC_FLAGS} -Wpadded")
   set (PEDANTIC_FLAGS ${DEFAULT_PEDANTIC_FLAGS} CACHE STRING "Compiler flags to enable pedantic warnings")

include/dohpgeom.h

@@ -57,7 +57,7 @@ INLINE void dGeomConvexComb_2_4(dReal x,dReal y,const dReal v[],const dInt p[],d
 * These are macros so we can use them with ints and entity handles.
 * 
 */
-#define dGeomMatch2(a0,a1,b0,b1) (((a0)==(b0) && (a1)==(b1)) || ((a0)==(b1) && (a1)==(b1)))
+#define dGeomMatch2(a0,a1,b0,b1) (((a0)==(b0) && (a1)==(b1)) || ((a0)==(b1) && (a1)==(b0)))
 #define dGeomMatchQuadLine(q,l)                                         \
   (dGeomMatch2((q)[0],(q)[1],(l)[0],(l)[1]) ? 0                         \
    : (dGeomMatch2((q)[1],(q)[2],(l)[0],(l)[1]) ? 1                      \

include/dohpmesh.h

@@ -115,12 +115,13 @@ EXTERN dErr dMeshGetNumEnts(dMesh,dMeshESH,dEntType,dEntTopology,dInt*);
 EXTERN dErr dMeshGetEnts(dMesh,dMeshESH,dEntType,dEntTopology,dMeshEH[],dInt,dInt*);
 EXTERN dErr dMeshGetAdjIndex(dMesh,const dMeshEH[],dInt,const dMeshEH[],dInt,dInt[],dInt*);
 
+EXTERN dErr dMeshGetTag(dMesh mesh,const char name[],dMeshTag *intag);
 EXTERN dErr dMeshTagDestroy(dMesh mesh,dMeshTag tag);
+EXTERN dErr dMeshTagCreate(dMesh mesh,const char[],dInt count,dDataType type,dMeshTag *intag);
 EXTERN dErr dMeshTagCreateTemp(dMesh mesh,const char[],dInt count,dDataType type,dMeshTag *intag);
 EXTERN dErr dMeshTagSetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,const void *data,dInt count,dDataType type);
 EXTERN dErr dMeshTagGetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,void *data,dInt count,dDataType type);
 EXTERN dErr dMeshGetTopo(dMesh,const dMeshEH[],dInt,dEntTopology[]);
-
 EXTERN dErr dMeshSetFromOptions(dMesh);
 EXTERN dErr dMeshTagBcast(dMesh mesh,dMeshTag tag);
 

include/dohptype.h

@@ -14,10 +14,13 @@ typedef PetscInt       dInt;
 typedef PetscReal      dReal;
 typedef PetscScalar    dScalar;
 typedef PetscTruth     dBool;
+typedef PetscTruth     dTruth;
 typedef PetscErrorCode dErr;
 typedef PetscObject    dObject;
 typedef PetscViewer    dViewer;
 
+typedef PetscMPIInt dMPIInt;
+
 /* #define dEntTopology      enum iMesh_EntityTopology */
 typedef int dMeshInt;
 typedef double dMeshReal;
@@ -29,7 +32,6 @@ typedef iBase_EntitySetHandle dMeshESH;
 typedef int dIInt;
 typedef double dIReal;
 typedef char dIByte;
-#define dFree7(a,b,c,d,e,f,g) PetscFree7((a),(b),(c),(d),(e),(f),(g))
 
 typedef enum { dDATA_INT, dDATA_REAL, dDATA_EH, dDATA_BYTE } dDataType;
 typedef enum { dTOPO_POINT, dTOPO_LINE, dTOPO_POLYGON, dTOPO_TRIANGLE,
@@ -48,7 +50,11 @@ typedef enum { dTYPE_VERTEX, dTYPE_EDGE, dTYPE_FACE, dTYPE_REGION, dTYPE_ALL } d
 #define dMemcpy(a,b,c) PetscMemcpy(a,b,c)
 #define dMemzero(a,b)  PetscMemzero(a,b)
 #define dValidHeader(a,b,c) PetscValidHeaderSpecific(a,b,c)
-#define dValidPointer(a,b) PetscValidPointer(a,b)
+
+#define dValidPointer(a,b)                                              \
+  { if (!(a)) dERROR(PETSC_ERR_ARG_NULL,"Null Pointer: Parameter # %d",(b)); \
+    if ((size_t)a & 3) dERROR(PETSC_ERR_ARG_BADPTR,"Invalid Pointer: Parameter # %d",(b));} \
+
 #define dMalloc(a,b) PetscMalloc(a,b)
 #define dNew(a,b) PetscNew(a,b)
 #define dNewLog(a,b,c) PetscNewLog(a,b,c)
@@ -97,10 +103,13 @@ typedef enum { dTYPE_VERTEX, dTYPE_EDGE, dTYPE_FACE, dTYPE_REGION, dTYPE_ALL } d
 #define dValidPointerSpecific7(p0,t0,a0,p1,t1,a1,p2,t2,a2,p3,t3,a3,p4,t4,a4,p5,t5,a5,p6,t6,a6) \
   {dValidPointerSpecific(p0,t0,a0); dValidPointerSpecific6(p1,t1,a1,p2,t2,a2,p3,t3,a3,p4,t4,a4,p5,t5,a5,p6,t6,a6);}
 
-#define dValidIntPointer(p,a) PetscValidIntPointer(p,a)
-#define dValidHandlePointer(p,a) dValidPointerNamedSpecific(p,void*,"void*",a)
-#define dValidScalarPointer(p,a) dValidPointerNamedSpecific(p,dScalar,"dScalar",a)
-#define dValidRealPointer(p,a) dValidPointerNamedSpecific(p,dReal,"dReal",a)
+#define dValidCharPointer(p,a) PetscValidCharPointer((p),(a))
+#define dValidIntPointer(p,a) PetscValidIntPointer((p),(a))
+#define dValidHandlePointer(p,a) dValidPointerNamedSpecific((p),void*,"void*",(a))
+#define dValidScalarPointer(p,a) dValidPointerNamedSpecific((p),dScalar,"dScalar",(a))
+#define dValidRealPointer(p,a) dValidPointerNamedSpecific((p),dReal,"dReal",(a))
+
+#define dStrlen(s,l) PetscStrlen((s),(l))
 
 #define dMax(a,b) PetscMax(a,b)
 #define dMin(a,b) PetscMin(a,b)
@@ -115,6 +124,10 @@ typedef enum { dTYPE_VERTEX, dTYPE_EDGE, dTYPE_FACE, dTYPE_REGION, dTYPE_ALL } d
 # define true PETSC_TRUE
 #endif
 
+/* stdbool has small (1 byte) bools, but ours are the same size as int to preserve alignment */
+#define dTRUE  PETSC_TRUE
+#define dFALSE PETSC_FALSE
+
 #define dMAX_PATH_LEN PETSC_MAX_PATH_LEN
 #define dNAME_LEN     256
 #define dSTR_LEN      256

sandbox/itaps/dohpblock.c

@@ -118,12 +118,11 @@ static dErr doMaterial(iMesh_Instance mesh)
     for (dInt j=0; j<r.s; j++) {
       dGeomVecMeanI(8,x.v+3*rvo.v[j],center);
       fx = sqrt(dGeomDotProd(center,center)); /* material 0 if inside the unit ball, else material 1 */
-      if (fx < 1.5) {
-        if (i == 0) ents[nents++] = r.v[j];
-      } else {
-        if (i == 1) ents[nents++] = r.v[j];
+      if (i == (fx < 1.0) ? 0 : 1) {
+        ents[nents] = r.v[j];
+        matnum[nents] = 1.0 * i;
+        nents++;
       }
-      matnum[j] = 1.0 * i;
     }
     iMesh_addEntArrToSet(mesh,ents,nents,&mat[i],&err);dICHK(mesh,err);
     iMesh_setArrData(mesh,ents,nents,matNumTag,(char*)matnum,nents*(int)sizeof(matnum[0]),&err);dICHK(mesh,err);
@@ -154,6 +153,31 @@ static dErr doGlobalNumber(iMesh_Instance mesh)
   dFunctionReturn(0);
 }
 
+static dErr doGlobalID(iMesh_Instance mesh)
+{
+  MeshListEH ents=MLZ;
+  MeshListInt type=MLZ;
+  int count[4] = {0,0,0,0};
+  int owned,offset,*number;
+  dMeshTag idTag;
+  dErr err;
+
+  dFunctionBegin;
+  iMesh_getEntities(mesh,0,iBase_ALL_TYPES,iMesh_ALL_TOPOLOGIES,MLREF(ents),&err);dICHK(mesh,err);
+  iMesh_getEntArrType(mesh,ents.v,ents.s,MLREF(type),&err);ICHKERRQ(mesh,err);
+  err = dMalloc(ents.s*sizeof(number[0]),&number);dCHK(err);
+  owned = ents.s; offset = 0;
+  for (int i=0; i<owned; i++) {
+    number[i] = count[type.v[i]]++;
+  }
+  iMesh_getTagHandle(mesh,"GLOBAL_ID",&idTag,&err,strlen("GLOBAL_ID"));dICHK(mesh,err);
+  iMesh_setIntArrData(mesh,ents.v,owned,idTag,number,owned,&err);dICHK(mesh,err);
+  err = dFree(number);dCHK(err);
+  MeshListFree(ents); MeshListFree(type);
+  dFunctionReturn(0);
+}
+
+
 static dErr createUniformTags(iMesh_Instance mesh)
 {
   dMeshTag itag,rtag;
@@ -183,8 +207,8 @@ static dErr createUniformTags(iMesh_Instance mesh)
 #define __FUNCT__ "main"
 int main(int argc, char *argv[])
 {
-  const char outopts[]="",pTagName[]="dohp_partition", pSetName[]="PARALLEL_PARTITION";
-  PetscTruth do_bdy = 0,do_material = 1,do_uniform = 1,do_global_number = 0;
+  const char outopts[]="",pTagName[]="OWNING_PART", pSetName[]="PARALLEL_PARTITION";
+  PetscTruth do_bdy = 0,do_material = 1,do_uniform = 1,do_global_number = 0,do_global_id = 1;
   PetscTruth do_partition = 1,do_faces = 1,do_edges = 1,do_orient = 0;
   char outfile[256] = "dblock.h5m";
   dInt verbose = 1;
@@ -209,7 +233,8 @@ int main(int argc, char *argv[])
     err = PetscOptionsTruth("-do_bdy","create boundary sets","none",do_bdy,&do_bdy,NULL);dCHK(err);
     err = PetscOptionsTruth("-do_material","create material sets","none",do_material,&do_material,NULL);dCHK(err);
     err = PetscOptionsTruth("-do_uniform","create uniform sets","none",do_uniform,&do_uniform,NULL);dCHK(err);
-    err = PetscOptionsTruth("-do_global_number","create global_number sets","none",do_global_number,&do_global_number,NULL);dCHK(err);
+    err = PetscOptionsTruth("-do_global_number","create global_number tags","none",do_global_number,&do_global_number,NULL);dCHK(err);
+    err = PetscOptionsTruth("-do_global_id","create GLOBAL_ID tags","none",do_global_id,&do_global_id,NULL);dCHK(err);
     err = PetscOptionsTruth("-do_partition","create partition sets","none",do_partition,&do_partition,NULL);dCHK(err);
     err = PetscOptionsTruth("-do_faces","create face entities","none",do_faces,&do_faces,NULL);dCHK(err);
     err = PetscOptionsTruth("-do_edges","create face entities","none",do_edges,&do_edges,NULL);dCHK(err);
@@ -270,6 +295,7 @@ int main(int argc, char *argv[])
   printf("region size %d, status size %d\n",r.s,s.s);
 
   if (do_global_number) {err = doGlobalNumber(mesh);dCHK(err);}
+  if (do_global_id) {err = doGlobalID(mesh);dCHK(err);}
 
   if (do_partition) {           /* Partition tags */
     /* Create partition. */
@@ -279,7 +305,7 @@ int main(int argc, char *argv[])
         for (k=0; k<p-1; k++) {
           I = i*M/(m-1); J = j*N/(n-1); K = k*P/(p-1);
           const int pind = (I*N+J)*P+K,eind = (i*(n-1)+j)*(p-1)+k;
-          part.v[eind] = pind + 1; /* Add 1 because MATERIAL_SET counts from 1 */
+          part.v[eind] = pind;
         }
       }
     }
@@ -316,7 +342,7 @@ int main(int argc, char *argv[])
           }
           printf("part[%d (%d,%d,%d)] has %d regions\n",(i*N+j)*P+k,i,j,k,(int)(entp-entbuf));
           iMesh_addEntArrToSet(mesh,entbuf,(int)(entp-entbuf),&partset,&err);dICHK(mesh,err);
-          iMesh_setEntSetIntData(mesh,partset,pTag,(i*N+j)*P+k+1,&err);dICHK(mesh,err);
+          iMesh_setEntSetIntData(mesh,partset,pTag,(i*N+j)*P+k,&err);dICHK(mesh,err);
         }
       }
     }

src/fs/impls/cont/cont.c

@@ -64,19 +64,25 @@ static dErr dFSContPropogateDegree(dFS fs)
   dMeshTag indexTag;
   dIInt ierr;
   dEntTopology topo;
-  dInt *eind,*deg,type,cnt,total,*aind,nents[4],estart[4];
+  dInt *eind,*deg,type,cnt,total,nvtx,*aind,nents[4],estart[4];
+  dMPIInt rank;
   dErr err;
 
   dFunctionBegin;
   dValidHeader(fs,dFS_COOKIE,1);
   mesh = fs->mesh;
   err = dMeshGetInstance(mesh,&mi);dCHK(err);
+  err = MPI_Comm_rank(((dObject)mesh)->comm,&rank);dCHK(err);
 
   /* Number all entities */
   err = dMeshTagCreateTemp(mesh,"index",1,dDATA_INT,&indexTag);dCHK(err);
-  err = dMeshGetNumEnts(mesh,fs->active,iBase_ALL_TYPES,iMesh_ALL_TOPOLOGIES,&total);dCHK(err);
+  err = dMeshGetNumEnts(mesh,fs->active,dTYPE_ALL,dTOPO_ALL,&total);dCHK(err);
+  err = dMeshGetNumEnts(mesh,fs->active,dTYPE_VERTEX,dTOPO_ALL,&nvtx);dCHK(err);
+  total -= nvtx;
   err = PetscMalloc2(total,dMeshEH,&ents,total,dInt,&eind);dCHK(err);
-  for (type=iBase_VERTEX,cnt=0; type<iBase_ALL_TYPES; type++) {
+  /* We don't propogate to vertices so we'll ignore them, but not preserve normal indexing */
+  nents[dTYPE_VERTEX] = 0; estart[dTYPE_VERTEX] = 0;
+  for (type=dTYPE_EDGE,cnt=0; type<dTYPE_ALL; type++) { /* Get out all non-vertex entities */
     err = dMeshGetEnts(mesh,fs->active,type,iMesh_ALL_TOPOLOGIES,ents+cnt,total-cnt,&nents[type]);dCHK(err);
     estart[type] = cnt;
     cnt += nents[type];
@@ -96,25 +102,28 @@ static dErr dFSContPropogateDegree(dFS fs)
   for (dInt i=0; i<total; i++) {
     eind[i] = i;
   }
-  err = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"nents=(%d %d %d %d)(%d) estart=(%d %d %d %d)\n",nents[0],nents[1],nents[2],nents[3],cnt,estart[0],estart[1],estart[2],estart[3]);
+  err = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] nents=(%d %d %d %d)(%d) estart=(%d %d %d %d)\n",rank,nents[0],nents[1],nents[2],nents[3],cnt,estart[0],estart[1],estart[2],estart[3]);
   err = dMeshTagSetData(mesh,indexTag,ents,total,eind,total,dDATA_INT);dCHK(err);
 
   err = dMalloc(3*total*sizeof(*deg),&deg);dCHK(err);
   err = dMeshTagGetData(mesh,fs->degreetag,ents,total,deg,3*total,dDATA_INT);dCHK(err); /* Get degree everywhere */
-  iMesh_getEntArrAdj(mi,ents,total,iBase_VERTEX,MLREF(conn),MLREF(connoff),&ierr);dICHK(mi,ierr); /* connectivity */
-  for (type=iBase_REGION; type>iBase_VERTEX; type--) { /* Propogate degrees downward to facets (and vertices) */
+  /* connectivity for all but vertices */
+  iMesh_getEntArrAdj(mi,ents+nents[dTYPE_VERTEX],total-nents[dTYPE_VERTEX],iBase_VERTEX,MLREF(conn),MLREF(connoff),&ierr);dICHK(mi,ierr);
+  for (type=iBase_REGION; type>iBase_EDGE; type--) { /* Propogate degrees downward to facets (not vertices) */
     switch (type) {
       case iBase_REGION: topo = iMesh_HEXAHEDRON; break;
       case iBase_FACE: topo = iMesh_QUADRILATERAL; break;
       case iBase_EDGE: topo = iMesh_LINE_SEGMENT; break;
       default: dERROR(1,"should not happen");
     }
+    /* Get adjacent entities */
     iMesh_getEntArrAdj(mi,ents+estart[type],nents[type],type-1,MLREF(adj),MLREF(adjoff),&err);dICHK(mi,err);
-    iMesh_getEntArrAdj(mi,adj.v,adj.s,iBase_VERTEX,MLREF(aconn),MLREF(aconnoff),&ierr);dICHK(mi,ierr); /* adjacent connectivity */
+    /* And adjacent connectivity */
+    iMesh_getEntArrAdj(mi,adj.v,adj.s,iBase_VERTEX,MLREF(aconn),MLREF(aconnoff),&ierr);dICHK(mi,ierr);
     err = dMalloc(adj.s*sizeof(*aind),&aind);dCHK(err);
     err = dMeshTagGetData(mesh,indexTag,adj.v,adj.s,aind,adj.s,dDATA_INT);dCHK(err); /* get indices of adj ents */
     for (dInt i=0; i<nents[type]; i++) {
-      err = dJacobiPropogateDown(fs->jacobi,topo,conn.v+connoff.v[estart[type]+i],deg+estart[type]+3*i,
+      err = dJacobiPropogateDown(fs->jacobi,topo,conn.v+connoff.v[estart[type]+i-nents[dTYPE_VERTEX]],deg+(estart[type]+i)*3,
                                  aconn.v+aconnoff.v[adjoff.v[i]],aind+adjoff.v[i],deg);dCHK(err);
     }
     MeshListFree(adj); MeshListFree(adjoff);

src/fs/interface/fs.c

@@ -110,11 +110,11 @@ dErr dFSDestroy(dFS fs)
   dFunctionReturn(0);
 }
 
-/** 
+/**
 * Builds a function space.  Enforcement of constraints is implementation dependent.
-* 
+*
 * @param fs the function space
-* 
+*
 * @return err
 */
 dErr dFSBuildSpace(dFS fs)
@@ -126,8 +126,8 @@ dErr dFSBuildSpace(dFS fs)
   dFunctionBegin;
   dValidHeader(fs,dFS_COOKIE,1);
   /**
-  * FIXME: 
-  * 
+  * FIXME:
+  *
   */
   err = dMeshGetInstance(mesh,&mi);dCHK(err);
   if (fs->ops->buildspace) {

src/fs/mesh/impls/pack/pack.c

@@ -62,7 +62,7 @@ typedef struct {
 */
 typedef struct {
   dInt nr,ns;
-  char *partition;
+  char partition[dSTR_LEN];
   PackData *owned;
   PackData *unowned;
 } Mesh_Pack;
@@ -599,23 +599,20 @@ static dErr dMeshDestroy_Pack(dMesh mesh)
   err = PackDataDestroy(pack->unowned);dCHK(err);
   err = PetscFree(pack->owned);dCHK(err);
   err = PetscFree(pack->unowned);dCHK(err);
-  err = dFree(mesh->data);dCHK(err);
+  err = dFree(mesh->data);dCHK(err); /* Like this so the pointer is zeroed */
   dFunctionReturn(0);
 }
 
 static dErr dMeshSetFromOptions_Pack(dMesh mesh)
 {
+  static const char defaultPartition[] = "PARALLEL_PARTITION";
   Mesh_Pack *pack = mesh->data;
-  char str[dSTR_LEN];
   dBool flg;
   dErr err;
 
   dFunctionBegin;
-  err = PetscOptionsString("-dmesh_partition","Name of partition tag","dMeshSetInFile",pack->partition,str,sizeof(str),&flg);dCHK(err);
-  if (flg) {
-    err = PetscStrfree(pack->partition);dCHK(err);
-    err = PetscStrallocpy(str,&pack->partition);dCHK(err);
-  }
+  if (!pack->partition[0]) {err = PetscStrcpy(pack->partition,defaultPartition);dCHK(err);}
+  err = PetscOptionsString("-dmesh_partition","Name of partition tag","dMeshSetInFile",defaultPartition,pack->partition,sizeof(pack->partition),&flg);dCHK(err);
   dFunctionReturn(0);
 }
 
@@ -630,7 +627,6 @@ dErr dMeshCreate_Pack(dMesh mesh)
   mesh->data = (void*)pack;
   err = PetscNew(PackData,&pack->owned);dCHK(err);
   err = PetscNew(PackData,&pack->unowned);dCHK(err);
-  err = PetscStrallocpy("PARALLEL_PARTITION",&pack->partition);dCHK(err);
 
   iMesh_newMesh("PARALLEL",&mesh->mi,&err,(int)strlen("PARALLEL"));dICHK(mesh->mi,err);