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