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: {
 13:   PetscFunctionBeginUser;
 14:   options->filename[0] = '\0';
 15:   options->interpolate = PETSC_FALSE;
 16:   options->meshNum     = 0;

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

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

 30:   PetscFunctionBegin;
 31:   dim = 3;
 32:   PetscCall(DMCreate(comm, dm));
 33:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Simple Hybrid Mesh"));
 34:   PetscCall(DMSetType(*dm, DMPLEX));
 35:   PetscCall(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:     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
 45:     if (interpolate) {
 46:       DM idm;

 48:       PetscCall(DMPlexInterpolate(*dm, &idm));
 49:       PetscCall(DMDestroy(dm));
 50:       *dm = idm;
 51:     }
 52:     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 53:   }
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 80:   dim = 3;
 81:   PetscCall(DMCreate(comm, dm));
 82:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Reverse Hybrid Mesh"));
 83:   PetscCall(DMSetType(*dm, DMPLEX));
 84:   PetscCall(DMSetDimension(*dm, dim));
 85:   {
 86:     /* Simple mesh with 2 hexes and 3 wedges */
 87:     PetscInt    numPoints[2]         = {16, 5};
 88:     PetscInt    coneSize[21]         = {8, 8, 6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 89:     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};
 90:     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};
 91:     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,
 92:                                         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};

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

 98:       PetscCall(DMPlexInterpolate(*dm, &idm));
 99:       PetscCall(DMDestroy(dm));
100:       *dm = idm;
101:     }
102:     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
103:   }
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

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

114:   PetscFunctionBegin;
115:   PetscCall(DMGetDimension(*dm, &dim));
116:   PetscCheck(dim == 3, PetscObjectComm((PetscObject)*dm), PETSC_ERR_SUP, "No support for dimension %" PetscInt_FMT, dim);
117:   PetscCall(DMPlexGetChart(*dm, &pStart, &pEnd));
118:   PetscCall(PetscMalloc1(pEnd - pStart, &ind));
119:   for (p = 0; p < pEnd - pStart; ++p) ind[p] = p;
120:   PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
121:   for (c = cStart; c < cEnd; ++c) {
122:     PetscInt coneSize;

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

132:     PetscCall(DMPlexGetConeSize(*dm, c, &coneSize));
133:     if (coneSize == 6) ind[c] = off[1]++;
134:     else ind[c] = off[0]++;
135:   }
136:   PetscCheck(off[0] == cEnd - Nhyb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of normal cells %" PetscInt_FMT " should be %" PetscInt_FMT, off[0], cEnd - Nhyb);
137:   PetscCheck(off[1] == cEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of hybrid cells %" PetscInt_FMT " should be %" PetscInt_FMT, off[1] - off[0], Nhyb);
138:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, ind, PETSC_OWN_POINTER, &perm));
139:   PetscCall(DMPlexPermute(*dm, perm, &pdm));
140:   PetscCall(ISDestroy(&perm));
141:   PetscCall(DMDestroy(dm));
142:   *dm = pdm;
143:   PetscFunctionReturn(PETSC_SUCCESS);
144: }

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

153:   PetscFunctionBegin;
154:   PetscCall(PetscStrlen(filename, &len));
155:   if (len) {
156:     PetscCall(DMPlexCreateFromFile(comm, filename, "ex34_plex", PETSC_FALSE, dm));
157:     PetscCall(OrderHybridMesh(dm));
158:     if (interpolate) {
159:       DM idm;

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

182: int main(int argc, char **argv)
183: {
184:   DM     dm;
185:   AppCtx user;

187:   PetscFunctionBeginUser;
188:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
189:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
190:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
191:   PetscCall(DMDestroy(&dm));
192:   PetscCall(PetscFinalize());
193:   return 0;
194: }

196: /*TEST

198:   test:
199:     suffix: 0
200:     args: -interpolate -dm_view ascii::ascii_info_detail

202:   # Test needs to be reworked
203:   test:
204:     TODO: broken
205:     suffix: 1
206:     args: -mesh_num 1 -interpolate -dm_view ascii::ascii_info_detail

208: TEST*/