Actual source code: ex16.c
1: static char help[] = "Tests for creation of submeshes\n\n";
3: #include <petscdmplex.h>
5: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
6: {
7: PetscFunctionBeginUser;
8: PetscCall(DMCreate(comm, dm));
9: PetscCall(DMSetType(*dm, DMPLEX));
10: PetscCall(DMSetFromOptions(*dm));
11: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: // Label half of the cells
16: static PetscErrorCode CreateHalfCellsLabel(DM dm, PetscBool lower, DMLabel *label)
17: {
18: PetscInt cStart, cEnd, cStartSub, cEndSub;
20: PetscFunctionBeginUser;
21: PetscCall(DMCreateLabel(dm, "cells"));
22: PetscCall(DMGetLabel(dm, "cells", label));
23: PetscCall(DMLabelClearStratum(*label, 1));
24: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
25: if (lower) {
26: cStartSub = cStart;
27: cEndSub = cEnd / 2;
28: } else {
29: cStartSub = cEnd / 2;
30: cEndSub = cEnd;
31: }
32: for (PetscInt c = cStartSub; c < cEndSub; ++c) PetscCall(DMLabelSetValue(*label, c, 1));
33: PetscCall(DMPlexLabelComplete(dm, *label));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: // Label everything on the right half of the domain
38: static PetscErrorCode CreateHalfDomainLabel(DM dm, PetscBool lower, PetscReal height, DMLabel *label)
39: {
40: PetscReal centroid[3];
41: PetscInt cStart, cEnd, cdim;
43: PetscFunctionBeginUser;
44: PetscCall(DMGetCoordinatesLocalSetUp(dm));
45: PetscCall(DMCreateLabel(dm, "cells"));
46: PetscCall(DMGetLabel(dm, "cells", label));
47: PetscCall(DMLabelClearStratum(*label, 1));
48: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
49: PetscCall(DMGetCoordinateDim(dm, &cdim));
50: for (PetscInt c = cStart; c < cEnd; ++c) {
51: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
52: if (height > 0.0 && PetscAbsReal(centroid[1] - height) > PETSC_SMALL) continue;
53: if (lower) {
54: if (centroid[0] < 0.5) PetscCall(DMLabelSetValue(*label, c, 1));
55: } else {
56: if (centroid[0] > 0.5) PetscCall(DMLabelSetValue(*label, c, 1));
57: }
58: }
59: PetscCall(DMPlexLabelComplete(dm, *label));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: // Create a line of faces at a given x value
64: static PetscErrorCode CreateLineLabel(DM dm, PetscReal x, DMLabel *label)
65: {
66: PetscReal centroid[3];
67: PetscInt fStart, fEnd;
69: PetscFunctionBeginUser;
70: PetscCall(DMGetCoordinatesLocalSetUp(dm));
71: PetscCall(DMCreateLabel(dm, "faces"));
72: PetscCall(DMGetLabel(dm, "faces", label));
73: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
74: for (PetscInt f = fStart; f < fEnd; ++f) {
75: PetscCall(DMPlexComputeCellGeometryFVM(dm, f, NULL, centroid, NULL));
76: if (PetscAbsReal(centroid[0] - x) < PETSC_SMALL) PetscCall(DMLabelSetValue(*label, f, 1));
77: }
78: PetscCall(DMPlexLabelComplete(dm, *label));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode CreateVolumeSubmesh(DM dm, PetscBool domain, PetscBool lower, PetscReal height, DM *subdm)
83: {
84: DMLabel label, map;
86: PetscFunctionBegin;
87: if (domain) PetscCall(CreateHalfDomainLabel(dm, lower, height, &label));
88: else PetscCall(CreateHalfCellsLabel(dm, lower, &label));
89: PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm));
90: PetscCall(PetscObjectSetName((PetscObject)*subdm, "Submesh"));
91: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subdm, "sub_"));
92: PetscCall(DMViewFromOptions(*subdm, NULL, "-dm_view"));
93: PetscCall(DMPlexGetSubpointMap(*subdm, &map));
94: PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view"));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode TestBoundaryField(DM dm)
99: {
100: DM subdm;
101: DMLabel label, map;
102: PetscFV fvm;
103: Vec gv;
104: PetscInt cdim;
106: PetscFunctionBeginUser;
107: PetscCall(CreateLineLabel(dm, 0.5, &label));
108: PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdm));
109: PetscCall(PetscObjectSetName((PetscObject)subdm, "Submesh"));
110: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subdm, "sub_"));
111: PetscCall(DMViewFromOptions(subdm, NULL, "-dm_view"));
112: PetscCall(DMPlexGetSubpointMap(subdm, &map));
113: PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view"));
115: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)subdm), &fvm));
116: PetscCall(PetscObjectSetName((PetscObject)fvm, "testField"));
117: PetscCall(PetscFVSetNumComponents(fvm, 1));
118: PetscCall(DMGetCoordinateDim(subdm, &cdim));
119: PetscCall(PetscFVSetSpatialDimension(fvm, cdim));
120: PetscCall(PetscFVSetFromOptions(fvm));
122: PetscCall(DMAddField(subdm, NULL, (PetscObject)fvm));
123: PetscCall(PetscFVDestroy(&fvm));
124: PetscCall(DMCreateDS(subdm));
126: PetscCall(DMCreateGlobalVector(subdm, &gv));
127: PetscCall(PetscObjectSetName((PetscObject)gv, "potential"));
128: PetscCall(VecSet(gv, 1.));
129: PetscCall(VecViewFromOptions(gv, NULL, "-vec_view"));
130: PetscCall(VecDestroy(&gv));
131: PetscCall(DMDestroy(&subdm));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: int main(int argc, char **argv)
136: {
137: DM dm, subdm;
138: PetscReal height = -1.0;
139: PetscBool volume = PETSC_TRUE, domain = PETSC_FALSE;
141: PetscFunctionBeginUser;
142: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
143: PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &height, NULL));
144: PetscCall(PetscOptionsGetBool(NULL, NULL, "-volume", &volume, NULL));
145: PetscCall(PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL));
147: PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
148: if (volume) {
149: PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_TRUE, height, &subdm));
150: PetscCall(DMSetFromOptions(subdm));
151: PetscCall(DMDestroy(&subdm));
152: PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_FALSE, height, &subdm));
153: PetscCall(DMSetFromOptions(subdm));
154: PetscCall(DMDestroy(&subdm));
155: } else {
156: PetscCall(TestBoundaryField(dm));
157: }
158: PetscCall(DMDestroy(&dm));
159: PetscCall(PetscFinalize());
160: return 0;
161: }
163: /*TEST
165: test:
166: suffix: 0
167: requires: triangle
168: args: -dm_coord_space 0 -sub_dm_plex_check_all \
169: -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view
171: test:
172: suffix: 0_vtk
173: requires: triangle
174: args: -dm_coord_space 0 -sub_dm_plex_check_all \
175: -dm_view vtk: -sub_dm_view vtk: -map_view
177: # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes
178: testset:
179: nsize: 3
180: requires: parmetis
181: args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \
182: -sub_dm_plex_check_all -dm_view -sub_dm_view
184: test:
185: suffix: 1
187: test:
188: suffix: 2
189: args: -domain
191: # This set tests that global numberings can be made when some strata are missing on a process
192: testset:
193: nsize: 3
194: requires: hdf5
195: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
196: -sub_dm_plex_check_all -sub_dm_view hdf5:subdm.h5
198: test:
199: suffix: 3
200: args: -domain -height 0.625
202: # This set tests that global numberings can be made when some strata are missing on a process
203: testset:
204: nsize: 3
205: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
206: -sub_dm_plex_check_all -sub_dm_view {{vtk:subdm.vtk: vtk:subdm.vtu :subdm.txt :subdm_d.txt:ascii_info_detail}}
208: test:
209: suffix: 3_vtk
210: args: -domain -height 0.625
212: # This test checks whether filter can extract a lower-dimensional manifold and output a field on it
213: testset:
214: requires: hdf5
215: args: -volume 0 -dm_plex_simplex 0 -sub_dm_view hdf5:subdm.h5 -vec_view hdf5:subdm.h5::append
217: test:
218: suffix: surface_2d
219: args: -dm_plex_box_faces 5,5
221: test:
222: suffix: surface_3d
223: args: -dm_plex_box_faces 5,5,5
225: # This test checks that if cells are present in the SF, the dm is marked as having an overlap
226: test:
227: suffix: surface_2d_2
228: nsize: 3
229: args: -dm_plex_box_faces 5,5 -domain -height 0.625
231: TEST*/