1 static char help[] = "Tests interpolation and output of hybrid meshes\n\n";
2
3 #include <petscdmplex.h>
4
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;
10
ProcessOptions(MPI_Comm comm,AppCtx * options)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;
17
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 }
25
CreateHybridMesh(MPI_Comm comm,PetscBool interpolate,DM * dm)26 static PetscErrorCode CreateHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm)
27 {
28 PetscInt dim;
29
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};
43
44 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
45 if (interpolate) {
46 DM idm;
47
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 }
56
57 /*
58 This is not a valid mesh. We need to either change to tensor quad prisms or regular triangular prisms.
59
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 */
CreateReverseHybridMesh(MPI_Comm comm,PetscBool interpolate,DM * dm)75 static PetscErrorCode CreateReverseHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm)
76 {
77 PetscInt dim;
78
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};
93
94 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
95 if (interpolate) {
96 DM idm;
97
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 }
106
OrderHybridMesh(DM * dm)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];
113
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;
123
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;
131
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 }
145
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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;
152
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;
160
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 }
181
main(int argc,char ** argv)182 int main(int argc, char **argv)
183 {
184 DM dm;
185 AppCtx user;
186
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 }
195
196 /*TEST
197
198 test:
199 suffix: 0
200 args: -interpolate -dm_view ascii::ascii_info_detail
201
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
207
208 TEST*/
209