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*/