1: static char help[] = "Read in a mesh and test whether it is valid\n\n";

  3:  #include <petscdmplex.h>
  4: #if defined(PETSC_HAVE_CGNS)
  5: #undef I /* Very old CGNS stupidly uses I as a variable, which fails when using complex. Curse you idiot package managers */
  6: #include <cgnslib.h>
  7: #endif
  8: #if defined(PETSC_HAVE_EXODUSII)
  9: #include <exodusII.h>
 10: #endif

 12: typedef struct {
 13:   PetscBool interpolate;                  /* Generate intermediate mesh elements */
 14:   char      filename[PETSC_MAX_PATH_LEN]; /* Mesh filename */
 15:   PetscInt  dim;
 16:   PetscErrorCode (**bcFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 17: } AppCtx;

 19: static PetscErrorCode zero(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 20: {
 21:   PetscInt i;
 22:   for (i = 0; i < dim; ++i) u[i] = 0.0;
 23:   return 0;
 24: }

 26: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 27: {

 31:   options->interpolate = PETSC_TRUE;
 32:   options->filename[0] = '\0';
 33:   options->dim         = 2;
 34:   options->bcFuncs     = NULL;

 36:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
 37:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex2.c", options->interpolate, &options->interpolate, NULL);
 38:   PetscOptionsString("-filename", "The mesh file", "ex2.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
 39:   PetscOptionsInt("-dim", "The dimension of problem used for non-file mesh", "ex2.c", options->dim, &options->dim, NULL);
 40:   PetscOptionsEnd();
 41:   return(0);
 42: }

 44: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 45: {
 46:   size_t         len;

 50:   PetscStrlen(user->filename, &len);
 51:   if (!len) {
 52:     DMLabel  label;
 53:     PetscInt id = 1;

 55:     DMPlexCreateBoxMesh(comm, user->dim, 2, user->interpolate, dm);
 56:     /* Mark boundary and set BC */
 57:     DMCreateLabel(*dm, "boundary");
 58:     DMGetLabel(*dm, "boundary", &label);
 59:     DMPlexMarkBoundaryFaces(*dm, label);
 60:     DMPlexLabelComplete(*dm, label);
 61:     PetscMalloc1(1, &user->bcFuncs);
 62:     user->bcFuncs[0] = zero;
 63:     DMAddBoundary(*dm, DM_BC_ESSENTIAL, "wall", "boundary", 0, 0, NULL, (void (*)(void)) user->bcFuncs[0], 1, &id, user);
 64:   } else {
 65:     DMPlexCreateFromFile(comm, user->filename, user->interpolate, dm);
 66:   }
 67:   PetscObjectSetName((PetscObject) *dm, "Mesh");
 68:   DMSetFromOptions(*dm);
 69:   DMViewFromOptions(*dm, NULL, "-dm_view");
 70:   return(0);
 71: }

 73: static PetscErrorCode CheckMeshTopology(DM dm)
 74: {
 75:   PetscInt       dim, coneSize, cStart;
 76:   PetscBool      isSimplex;

 80:   DMGetDimension(dm, &dim);
 81:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
 82:   DMPlexGetConeSize(dm, cStart, &coneSize);
 83:   isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE;
 84:   DMPlexCheckSymmetry(dm);
 85:   DMPlexCheckSkeleton(dm, isSimplex, 0);
 86:   DMPlexCheckFaces(dm, isSimplex, 0);
 87:   return(0);
 88: }

 90: static PetscErrorCode CheckMeshGeometry(DM dm)
 91: {
 92:   PetscInt       dim, coneSize, cStart, cEnd, c;
 93:   PetscReal     *v0, *J, *invJ, detJ;

 97:   DMGetDimension(dm, &dim);
 98:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 99:   DMPlexGetConeSize(dm, cStart, &coneSize);
100:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
101:   for (c = cStart; c < cEnd; ++c) {
102:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
103:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for cell %D", (double)detJ, c);
104:   }
105:   PetscFree3(v0,J,invJ);
106:   return(0);
107: }

109: int main(int argc, char **argv)
110: {
111:   DM             dm;
112:   AppCtx         user;

115:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
116:   ProcessOptions(PETSC_COMM_WORLD, &user);
117:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
118:   CheckMeshTopology(dm);
119:   CheckMeshGeometry(dm);
120:   DMDestroy(&dm);
121:   PetscFree(user.bcFuncs);
122:   PetscFinalize();
123:   return ierr;
124: }

126: /*TEST

128:   # CGNS meshes 0-1
129:   test:
130:     suffix: 0
131:     requires: cgns
132:     TODO: broken
133:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/tut21.cgns -interpolate 1
134:   test:
135:     suffix: 1
136:     requires: cgns
137:     TODO: broken
138:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/grid_c.cgns -interpolate 1
139:   # Gmsh meshes 2-4
140:   test:
141:     suffix: 2
142:     requires: double
143:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -interpolate 1
144:   test:
145:     suffix: 3
146:     requires: double
147:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square.msh -interpolate 1
148:   test:
149:     suffix: 4
150:     requires: double
151:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_bin.msh -interpolate 1
152:   # Exodus meshes 5-9
153:   test:
154:     suffix: 5
155:     requires: exodusii
156:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -interpolate 1
157:   test:
158:     suffix: 6
159:     requires: exodusii
160:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -interpolate 1
161:   test:
162:     suffix: 7
163:     requires: exodusii
164:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/squaremotor-30.exo -interpolate 1
165:   test:
166:     suffix: 8
167:     requires: exodusii
168:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -interpolate 1
169:   test:
170:     suffix: 9
171:     requires: exodusii
172:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/simpleblock-100.exo -interpolate 1

174: TEST*/