Actual source code: ex57.c
1: static char help[] = "Tests for ephemeral meshes.\n";
3: #include <petscdmplex.h>
4: #include <petscdmplextransform.h>
5: #include <petscdmlabelephemeral.h>
7: /*
8: Use
10: -ref_dm_view -eph_dm_view
12: to view the concrete and ephemeral meshes from the first transformation, and
14: -ref_dm_sec_view -eph_dm_sec_view
16: for the second.
17: */
19: // Should remove when I have an API for everything
20: #include <petsc/private/dmplextransformimpl.h>
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;
28: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
29: {
30: PetscInt cells[1024], Nc = 1024;
31: PetscBool flg;
33: PetscFunctionBeginUser;
34: options->active = NULL;
35: options->second = PETSC_FALSE;
36: options->concrete = PETSC_FALSE;
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: }
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: }
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));
72: PetscCall(DMSetApplicationContext(dm, *tr));
73: PetscCall(PetscObjectViewFromOptions((PetscObject)*tr, NULL, "-dm_plex_transform_view"));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
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: }
85: static PetscErrorCode CreateConcreteMesh(DMPlexTransform tr, DM *rdm)
86: {
87: DM cdm, codm, rcodm;
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: }
110: // dmA is concrete and dmB is ephemeral
111: static PetscErrorCode CompareMeshes(DM dmA, DM dmB, DM dm)
112: {
113: PetscInt dim, dimB, pStart, pEnd, pStartB, pEndB;
114: MPI_Comm comm;
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;
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: }
142: int main(int argc, char *argv[])
143: {
144: DM dm, tdm, rdm;
145: DMPlexTransform tr;
146: AppCtx user;
147: MPI_Comm comm;
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;
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: }
185: /*TEST
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
193: test:
194: suffix: quad
195: args: -dm_plex_simplex 0
197: # Here I am checking that the 'marker' label is correct for the ephemeral mesh
198: test:
199: suffix: tet
200: requires: ctetgen
201: args: -dm_plex_dim 3 -ref_dm_view -eph_dm_view -eph_dm_plex_view_labels marker
203: test:
204: suffix: hex
205: args: -dm_plex_dim 3 -dm_plex_simplex 0
207: # Tests for filter patches
208: testset:
209: args: -first_dm_plex_transform_type transform_filter -ref_dm_view
211: test:
212: suffix: tri_patch
213: requires: triangle
214: args: -cells 0,1,2,4
216: test:
217: suffix: quad_patch
218: args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -cells 0,1,3,4
220: # Tests for refined filter patches
221: testset:
222: args: -first_dm_plex_transform_type transform_filter -ref_dm_view -eph_dm_view -second
224: test:
225: suffix: tri_patch_ref
226: requires: triangle
227: args: -cells 0,1,2,4
229: test:
230: suffix: tri_patch_ref_concrete
231: requires: triangle
232: args: -cells 0,1,2,4 -concrete -first_dm_plex_transform_view ::ascii_info_detail
234: # Tests for boundary layer refinement
235: test:
236: suffix: quad_bl
237: args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_extrude 2 -cells 0,2,4,6,8 \
238: -first_dm_plex_transform_type refine_boundary_layer -first_dm_plex_transform_bl_splits 4 \
239: -ref_dm_view
241: TEST*/