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