Actual source code: ex34.c

  1: static char help[] = "Tests interpolation and output of hybrid meshes\n\n";

  3: #include <petscdmplex.h>

  5: typedef struct {
  6:   char      filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
  7:   PetscBool interpolate;                  /* Interpolate the mesh */
  8:   PetscInt  meshNum;                      /* Which mesh we should construct */
  9: } AppCtx;

 11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 12: {
 14:   options->filename[0] = '\0';
 15:   options->interpolate = PETSC_FALSE;
 16:   options->meshNum     = 0;

 18:   PetscOptionsBegin(comm, "", "Hybrid Output Test Options", "DMPLEX");
 19:   PetscOptionsString("-filename", "The mesh file", "ex8.c", options->filename, options->filename, sizeof(options->filename), NULL);
 20:   PetscOptionsBool("-interpolate", "Interpolate the mesh", "ex8.c", options->interpolate, &options->interpolate, NULL);
 21:   PetscOptionsBoundedInt("-mesh_num", "The mesh we should construct", "ex8.c", options->meshNum, &options->meshNum, NULL, 0);
 22:   PetscOptionsEnd();

 24:   return 0;
 25: }

 27: static PetscErrorCode CreateHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm)
 28: {
 29:   PetscInt dim;

 31:   dim = 3;
 32:   DMCreate(comm, dm);
 33:   PetscObjectSetName((PetscObject)*dm, "Simple Hybrid Mesh");
 34:   DMSetType(*dm, DMPLEX);
 35:   DMSetDimension(*dm, dim);
 36:   {
 37:     /* Simple mesh with 2 tets and 1 wedge */
 38:     PetscInt    numPoints[2]         = {8, 3};
 39:     PetscInt    coneSize[11]         = {4, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0};
 40:     PetscInt    cones[14]            = {4, 5, 6, 3, 7, 9, 8, 10, 4, 5, 6, 7, 8, 9};
 41:     PetscInt    coneOrientations[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 42:     PetscScalar vertexCoords[48]     = {-1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 0.0};

 44:     DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 45:     if (interpolate) {
 46:       DM idm;

 48:       DMPlexInterpolate(*dm, &idm);
 49:       DMDestroy(dm);
 50:       *dm = idm;
 51:     }
 52:     DMViewFromOptions(*dm, NULL, "-dm_view");
 53:   }
 54:   return 0;
 55: }

 57: /*
 58:    This is not a valid mesh. We need to either change to tensor quad prisms or regular triangular prisms.

 60:            10-------16--------20
 61:            /|        |
 62:           / |        |
 63:          /  |        |
 64:         9---|---15   |
 65:        /|   7    |  13--------18
 66:       / |  /     |  /    ____/
 67:      /  | /      | /____/
 68:     8   |/  14---|//---19
 69:     |   6    |  12
 70:     |  /     |  / \
 71:     | /      | /   \__
 72:     |/       |/       \
 73:     5--------11--------17
 74: */
 75: static PetscErrorCode CreateReverseHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm)
 76: {
 77:   PetscInt dim;

 79:   dim = 3;
 80:   DMCreate(comm, dm);
 81:   PetscObjectSetName((PetscObject)*dm, "Reverse Hybrid Mesh");
 82:   DMSetType(*dm, DMPLEX);
 83:   DMSetDimension(*dm, dim);
 84:   {
 85:     /* Simple mesh with 2 hexes and 3 wedges */
 86:     PetscInt    numPoints[2]         = {16, 5};
 87:     PetscInt    coneSize[21]         = {8, 8, 6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 88:     PetscInt    cones[34]            = {5, 6, 12, 11, 8, 14, 15, 9, 6, 7, 13, 12, 9, 15, 16, 10, 11, 17, 12, 14, 19, 15, 12, 18, 13, 15, 20, 16, 12, 17, 18, 15, 19, 20};
 89:     PetscInt    coneOrientations[34] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 90:     PetscScalar vertexCoords[48]     = {-1.0, -1.0, 0.0, -1.0, 0.0,  0.0, -1.0, 1.0, 0.0, -1.0, -1.0, 1.0, -1.0, 0.0,  1.0, -1.0, 1.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0,
 91:                                         0.0,  1.0,  0.0, 0.0,  -1.0, 1.0, 0.0,  0.0, 1.0, 0.0,  1.0,  1.0, 1.0,  -1.0, 0.0, 1.0,  1.0, 0.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0};

 93:     DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 94:     if (interpolate) {
 95:       DM idm;

 97:       DMPlexInterpolate(*dm, &idm);
 98:       DMDestroy(dm);
 99:       *dm = idm;
100:     }
101:     DMViewFromOptions(*dm, NULL, "-dm_view");
102:   }
103:   return 0;
104: }

106: static PetscErrorCode OrderHybridMesh(DM *dm)
107: {
108:   DM        pdm;
109:   IS        perm;
110:   PetscInt *ind;
111:   PetscInt  dim, pStart, pEnd, p, cStart, cEnd, c, Nhyb = 0, off[2];

113:   DMGetDimension(*dm, &dim);
115:   DMPlexGetChart(*dm, &pStart, &pEnd);
116:   PetscMalloc1(pEnd - pStart, &ind);
117:   for (p = 0; p < pEnd - pStart; ++p) ind[p] = p;
118:   DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
119:   for (c = cStart; c < cEnd; ++c) {
120:     PetscInt coneSize;

122:     DMPlexGetConeSize(*dm, c, &coneSize);
123:     if (coneSize == 6) ++Nhyb;
124:   }
125:   off[0] = 0;
126:   off[1] = cEnd - Nhyb;
127:   for (c = cStart; c < cEnd; ++c) {
128:     PetscInt coneSize;

130:     DMPlexGetConeSize(*dm, c, &coneSize);
131:     if (coneSize == 6) ind[c] = off[1]++;
132:     else ind[c] = off[0]++;
133:   }
136:   ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, ind, PETSC_OWN_POINTER, &perm);
137:   DMPlexPermute(*dm, perm, &pdm);
138:   ISDestroy(&perm);
139:   DMDestroy(dm);
140:   *dm = pdm;
141:   return 0;
142: }

144: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
145: {
146:   const char *filename    = user->filename;
147:   PetscBool   interpolate = user->interpolate;
148:   PetscInt    meshNum     = user->meshNum;
149:   size_t      len;

151:   PetscStrlen(filename, &len);
152:   if (len) {
153:     DMPlexCreateFromFile(comm, filename, "ex34_plex", PETSC_FALSE, dm);
154:     OrderHybridMesh(dm);
155:     if (interpolate) {
156:       DM idm;

158:       DMPlexInterpolate(*dm, &idm);
159:       DMDestroy(dm);
160:       *dm = idm;
161:     }
162:     PetscObjectSetName((PetscObject)*dm, "Input Mesh");
163:     DMViewFromOptions(*dm, NULL, "-dm_view");
164:   } else {
165:     switch (meshNum) {
166:     case 0:
167:       CreateHybridMesh(comm, interpolate, dm);
168:       break;
169:     case 1:
170:       CreateReverseHybridMesh(comm, interpolate, dm);
171:       break;
172:     default:
173:       SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown mesh number %" PetscInt_FMT, user->meshNum);
174:     }
175:   }
176:   return 0;
177: }

179: int main(int argc, char **argv)
180: {
181:   DM     dm;
182:   AppCtx user;

185:   PetscInitialize(&argc, &argv, NULL, help);
186:   ProcessOptions(PETSC_COMM_WORLD, &user);
187:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
188:   DMDestroy(&dm);
189:   PetscFinalize();
190:   return 0;
191: }

193: /*TEST

195:   test:
196:     suffix: 0
197:     args: -interpolate -dm_view ascii::ascii_info_detail

199:   # Test needs to be reworked
200:   test:
201:     requires: BROKEN
202:     suffix: 1
203:     args: -mesh_num 1 -interpolate -dm_view ascii::ascii_info_detail

205: TEST*/