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();

 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

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

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

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

 49:       PetscCall(DMPlexInterpolate(*dm, &idm));
 50:       PetscCall(DMDestroy(dm));
 51:       *dm = idm;
 52:     }
 53:     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 54:   }
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

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

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

 80:   PetscFunctionBegin;
 81:   dim = 3;
 82:   PetscCall(DMCreate(comm, dm));
 83:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Reverse Hybrid Mesh"));
 84:   PetscCall(DMSetType(*dm, DMPLEX));
 85:   PetscCall(DMSetDimension(*dm, dim));
 86:   {
 87:     /* Simple mesh with 2 hexes and 3 wedges */
 88:     PetscInt    numPoints[2]         = {16, 5};
 89:     PetscInt    coneSize[21]         = {8, 8, 6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 90:     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};
 91:     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};
 92:     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,
 93:                                         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};

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

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

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

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

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

133:     PetscCall(DMPlexGetConeSize(*dm, c, &coneSize));
134:     if (coneSize == 6) ind[c] = off[1]++;
135:     else ind[c] = off[0]++;
136:   }
137:   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);
138:   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);
139:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, ind, PETSC_OWN_POINTER, &perm));
140:   PetscCall(DMPlexPermute(*dm, perm, &pdm));
141:   PetscCall(ISDestroy(&perm));
142:   PetscCall(DMDestroy(dm));
143:   *dm = pdm;
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

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

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

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

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

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

197: /*TEST

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

203:   # Test needs to be reworked
204:   test:
205:     requires: BROKEN
206:     suffix: 1
207:     args: -mesh_num 1 -interpolate -dm_view ascii::ascii_info_detail

209: TEST*/