1 static char help[] = "Tests for ephemeral meshes.\n";
2
3 #include <petscdmplex.h>
4 #include <petscdmplextransform.h>
5 #include <petscdmlabelephemeral.h>
6
7 /*
8 Use
9
10 -ref_dm_view -eph_dm_view
11
12 to view the concrete and ephemeral meshes from the first transformation, and
13
14 -ref_dm_sec_view -eph_dm_sec_view
15
16 for the second.
17 */
18
19 // Should remove when I have an API for everything
20 #include <petsc/private/dmplextransformimpl.h>
21
22 typedef struct {
23 DMLabel active; /* Label for transform */
24 PetscBool second; /* Flag to execute a second transformation */
25 PetscBool concrete; /* Flag to use the concrete mesh for the second transformation */
26 } AppCtx;
27
ProcessOptions(MPI_Comm comm,AppCtx * options)28 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
29 {
30 PetscInt cells[1024], Nc = 1024;
31 PetscBool flg;
32
33 PetscFunctionBeginUser;
34 options->active = NULL;
35 options->second = PETSC_FALSE;
36 options->concrete = PETSC_FALSE;
37
38 PetscOptionsBegin(comm, "", "Ephemeral Meshing Options", "DMPLEX");
39 PetscCall(PetscOptionsIntArray("-cells", "Cells to mark for transformation", "ex57.c", cells, &Nc, &flg));
40 if (flg) {
41 PetscCall(DMLabelCreate(comm, "active", &options->active));
42 for (PetscInt c = 0; c < Nc; ++c) PetscCall(DMLabelSetValue(options->active, cells[c], DM_ADAPT_REFINE));
43 }
44 PetscCall(PetscOptionsBool("-second", "Use a second transformation", "ex57.c", options->second, &options->second, &flg));
45 PetscCall(PetscOptionsBool("-concrete", "Use concrete mesh for the second transformation", "ex57.c", options->concrete, &options->concrete, &flg));
46 PetscOptionsEnd();
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
CreateMesh(MPI_Comm comm,DM * dm)50 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
51 {
52 PetscFunctionBeginUser;
53 PetscCall(DMCreate(comm, dm));
54 PetscCall(DMSetType(*dm, DMPLEX));
55 PetscCall(DMSetFromOptions(*dm));
56 PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
57 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
CreateTransform(DM dm,DMLabel active,const char prefix[],DMPlexTransform * tr)61 static PetscErrorCode CreateTransform(DM dm, DMLabel active, const char prefix[], DMPlexTransform *tr)
62 {
63 PetscFunctionBeginUser;
64 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), tr));
65 PetscCall(PetscObjectSetName((PetscObject)*tr, "Transform"));
66 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*tr, prefix));
67 PetscCall(DMPlexTransformSetDM(*tr, dm));
68 PetscCall(DMPlexTransformSetActive(*tr, active));
69 PetscCall(DMPlexTransformSetFromOptions(*tr));
70 PetscCall(DMPlexTransformSetUp(*tr));
71
72 PetscCall(DMSetApplicationContext(dm, *tr));
73 PetscCall(PetscObjectViewFromOptions((PetscObject)*tr, NULL, "-dm_plex_transform_view"));
74 PetscFunctionReturn(PETSC_SUCCESS);
75 }
76
CreateEphemeralMesh(DMPlexTransform tr,DM * tdm)77 static PetscErrorCode CreateEphemeralMesh(DMPlexTransform tr, DM *tdm)
78 {
79 PetscFunctionBegin;
80 PetscCall(DMPlexCreateEphemeral(tr, "eph_", tdm));
81 PetscCall(PetscObjectSetName((PetscObject)*tdm, "Ephemeral Mesh"));
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
CreateConcreteMesh(DMPlexTransform tr,DM * rdm)85 static PetscErrorCode CreateConcreteMesh(DMPlexTransform tr, DM *rdm)
86 {
87 DM cdm, codm, rcodm;
88
89 PetscFunctionBegin;
90 PetscCall(DMPlexTransformGetDM(tr, &cdm));
91 PetscCall(DMPlexTransformApply(tr, cdm, rdm));
92 PetscCall(DMSetCoarsenLevel(*rdm, cdm->leveldown));
93 PetscCall(DMSetRefineLevel(*rdm, cdm->levelup + 1));
94 PetscCall(DMCopyDisc(cdm, *rdm));
95 PetscCall(DMGetCoordinateDM(cdm, &codm));
96 PetscCall(DMGetCoordinateDM(*rdm, &rcodm));
97 PetscCall(DMCopyDisc(codm, rcodm));
98 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
99 PetscCall(DMSetCoarseDM(*rdm, cdm));
100 PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
101 if (rdm) {
102 ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)cdm->data)->printFEM;
103 ((DM_Plex *)(*rdm)->data)->printL2 = ((DM_Plex *)cdm->data)->printL2;
104 }
105 PetscCall(PetscObjectSetName((PetscObject)*rdm, "Concrete Mesh"));
106 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*rdm, "ref_"));
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
110 // dmA is concrete and dmB is ephemeral
CompareMeshes(DM dmA,DM dmB,DM dm)111 static PetscErrorCode CompareMeshes(DM dmA, DM dmB, DM dm)
112 {
113 PetscInt dim, dimB, pStart, pEnd, pStartB, pEndB;
114 MPI_Comm comm;
115
116 PetscFunctionBegin;
117 PetscCall(PetscObjectGetComm((PetscObject)dmA, &comm));
118 PetscCall(DMGetDimension(dmA, &dim));
119 PetscCall(DMGetDimension(dmB, &dimB));
120 PetscCheck(dim == dimB, comm, PETSC_ERR_ARG_INCOMP, "Dimension from dmA %" PetscInt_FMT " != %" PetscInt_FMT " from dmB", dim, dimB);
121 PetscCall(DMPlexGetChart(dmA, &pStart, &pEnd));
122 PetscCall(DMPlexGetChart(dmB, &pStartB, &pEndB));
123 PetscCheck(pStart == pStartB && pEnd == pEndB, comm, PETSC_ERR_ARG_INCOMP, "Chart from dmA (%" PetscInt_FMT ", %" PetscInt_FMT ") does not match (%" PetscInt_FMT ", %" PetscInt_FMT ") for dmB", pStart, pEnd, pStartB, pEndB);
124 for (PetscInt p = pStart; p < pEnd; ++p) {
125 const PetscInt *cone, *ornt, *coneB, *orntB;
126 PetscInt coneSize, coneSizeB;
127
128 PetscCall(DMPlexGetConeSize(dmA, p, &coneSize));
129 PetscCall(DMPlexGetConeSize(dmB, p, &coneSizeB));
130 PetscCheck(coneSize == coneSizeB, comm, PETSC_ERR_ARG_INCOMP, "Cone size for %" PetscInt_FMT " from dmA %" PetscInt_FMT " does not match %" PetscInt_FMT " for dmB", p, coneSize, coneSizeB);
131 PetscCall(DMPlexGetOrientedCone(dmA, p, &cone, &ornt));
132 PetscCall(DMPlexGetOrientedCone(dmB, p, &coneB, &orntB));
133 for (PetscInt c = 0; c < coneSize; ++c) {
134 PetscCheck(cone[c] == coneB[c] && ornt[c] == orntB[c], comm, PETSC_ERR_ARG_INCOMP, "Cone point %" PetscInt_FMT " for point %" PetscInt_FMT " from dmA (%" PetscInt_FMT ", %" PetscInt_FMT ") does not match (%" PetscInt_FMT ", %" PetscInt_FMT ") for dmB", c, p, cone[c], ornt[c], coneB[c], orntB[c]);
135 }
136 PetscCall(DMPlexRestoreOrientedCone(dmA, p, &cone, &ornt));
137 PetscCall(DMPlexRestoreOrientedCone(dmB, p, &coneB, &orntB));
138 }
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141
main(int argc,char * argv[])142 int main(int argc, char *argv[])
143 {
144 DM dm, tdm, rdm;
145 DMPlexTransform tr;
146 AppCtx user;
147 MPI_Comm comm;
148
149 PetscFunctionBeginUser;
150 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
151 comm = PETSC_COMM_WORLD;
152 PetscCall(ProcessOptions(comm, &user));
153 PetscCall(CreateMesh(comm, &dm));
154 PetscCall(CreateTransform(dm, user.active, "first_", &tr));
155 PetscCall(CreateEphemeralMesh(tr, &tdm));
156 PetscCall(CreateConcreteMesh(tr, &rdm));
157 if (user.second) {
158 DMPlexTransform tr2;
159 DM tdm2, rdm2;
160
161 PetscCall(DMViewFromOptions(rdm, NULL, "-dm_sec_view"));
162 PetscCall(DMViewFromOptions(tdm, NULL, "-dm_sec_view"));
163 if (user.concrete) PetscCall(CreateTransform(rdm, user.active, "second_", &tr2));
164 else PetscCall(CreateTransform(tdm, user.active, "second_", &tr2));
165 PetscCall(CreateEphemeralMesh(tr2, &tdm2));
166 PetscCall(CreateConcreteMesh(tr2, &rdm2));
167 PetscCall(DMDestroy(&tdm));
168 PetscCall(DMDestroy(&rdm));
169 PetscCall(DMPlexTransformDestroy(&tr2));
170 tdm = tdm2;
171 rdm = rdm2;
172 }
173 PetscCall(DMViewFromOptions(tdm, NULL, "-dm_view"));
174 PetscCall(DMViewFromOptions(rdm, NULL, "-dm_view"));
175 PetscCall(CompareMeshes(rdm, tdm, dm));
176 PetscCall(DMPlexTransformDestroy(&tr));
177 PetscCall(DMDestroy(&tdm));
178 PetscCall(DMDestroy(&rdm));
179 PetscCall(DMDestroy(&dm));
180 PetscCall(DMLabelDestroy(&user.active));
181 PetscCall(PetscFinalize());
182 return 0;
183 }
184
185 /*TEST
186
187 # Tests for regular refinement of whole meshes
188 test:
189 suffix: tri
190 requires: triangle
191 args: -first_dm_plex_transform_view ::ascii_info_detail
192
193 test:
194 suffix: quad
195 args: -dm_plex_simplex 0
196 output_file: output/empty.out
197
198 # Here I am checking that the 'marker' label is correct for the ephemeral mesh
199 test:
200 suffix: tet
201 requires: ctetgen
202 args: -dm_plex_dim 3 -ref_dm_view -eph_dm_view -eph_dm_plex_view_labels marker
203
204 test:
205 suffix: hex
206 args: -dm_plex_dim 3 -dm_plex_simplex 0
207 output_file: output/empty.out
208
209 # Tests for filter patches
210 testset:
211 args: -first_dm_plex_transform_type transform_filter -ref_dm_view
212
213 test:
214 suffix: tri_patch
215 requires: triangle
216 args: -cells 0,1,2,4
217
218 test:
219 suffix: quad_patch
220 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -cells 0,1,3,4
221
222 # Tests for refined filter patches
223 testset:
224 args: -first_dm_plex_transform_type transform_filter -ref_dm_view -eph_dm_view -second
225
226 test:
227 suffix: tri_patch_ref
228 requires: triangle
229 args: -cells 0,1,2,4
230
231 test:
232 suffix: tri_patch_ref_concrete
233 requires: triangle
234 args: -cells 0,1,2,4 -concrete -first_dm_plex_transform_view ::ascii_info_detail
235
236 # Tests for boundary layer refinement
237 test:
238 suffix: quad_bl
239 args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_extrude 2 -cells 0,2,4,6,8 \
240 -first_dm_plex_transform_type refine_boundary_layer -first_dm_plex_transform_bl_splits 4 \
241 -ref_dm_view
242
243 # Tests for extrusion
244 test:
245 suffix: sphere_extruded
246 args: -dm_plex_shape sphere \
247 -first_dm_plex_transform_type extrude \
248 -first_dm_plex_transform_extrude_layers 3 \
249 -first_dm_plex_transform_extrude_use_tensor 0
250 output_file: output/empty.out
251
252 TEST*/
253