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: {
8: DMCreate(comm, dm);
9: DMSetType(*dm, DMPLEX);
10: DMSetFromOptions(*dm);
11: DMViewFromOptions(*dm, NULL, "-dm_view");
12: return 0;
13: }
15: // Label half of the cells
16: static PetscErrorCode CreateHalfCellsLabel(DM dm, PetscBool lower, DMLabel *label)
17: {
18: PetscInt cStart, cEnd, cStartSub, cEndSub;
21: DMCreateLabel(dm, "cells");
22: DMGetLabel(dm, "cells", label);
23: DMLabelClearStratum(*label, 1);
24: 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) DMLabelSetValue(*label, c, 1);
33: DMPlexLabelComplete(dm, *label);
34: return 0;
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;
44: DMGetCoordinatesLocalSetUp(dm);
45: DMCreateLabel(dm, "cells");
46: DMGetLabel(dm, "cells", label);
47: DMLabelClearStratum(*label, 1);
48: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
49: DMGetCoordinateDim(dm, &cdim);
50: for (PetscInt c = cStart; c < cEnd; ++c) {
51: 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) DMLabelSetValue(*label, c, 1);
55: } else {
56: if (centroid[0] > 0.5) DMLabelSetValue(*label, c, 1);
57: }
58: }
59: DMPlexLabelComplete(dm, *label);
60: return 0;
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;
70: DMGetCoordinatesLocalSetUp(dm);
71: DMCreateLabel(dm, "faces");
72: DMGetLabel(dm, "faces", label);
73: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
74: for (PetscInt f = fStart; f < fEnd; ++f) {
75: DMPlexComputeCellGeometryFVM(dm, f, NULL, centroid, NULL);
76: if (PetscAbsReal(centroid[0] - x) < PETSC_SMALL) DMLabelSetValue(*label, f, 1);
77: }
78: DMPlexLabelComplete(dm, *label);
79: return 0;
80: }
82: static PetscErrorCode CreateVolumeSubmesh(DM dm, PetscBool domain, PetscBool lower, PetscReal height, DM *subdm)
83: {
84: DMLabel label, map;
86: if (domain) CreateHalfDomainLabel(dm, lower, height, &label);
87: else CreateHalfCellsLabel(dm, lower, &label);
88: DMPlexFilter(dm, label, 1, subdm);
89: PetscObjectSetName((PetscObject)*subdm, "Submesh");
90: PetscObjectSetOptionsPrefix((PetscObject)*subdm, "sub_");
91: DMViewFromOptions(*subdm, NULL, "-dm_view");
92: DMPlexGetSubpointMap(*subdm, &map);
93: PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view");
94: return 0;
95: }
97: static PetscErrorCode TestBoundaryField(DM dm)
98: {
99: DM subdm;
100: DMLabel label, map;
101: PetscFV fvm;
102: Vec gv;
103: PetscInt cdim;
106: CreateLineLabel(dm, 0.5, &label);
107: DMPlexFilter(dm, label, 1, &subdm);
108: PetscObjectSetName((PetscObject)subdm, "Submesh");
109: PetscObjectSetOptionsPrefix((PetscObject)subdm, "sub_");
110: DMViewFromOptions(subdm, NULL, "-dm_view");
111: DMPlexGetSubpointMap(subdm, &map);
112: PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view");
114: PetscFVCreate(PetscObjectComm((PetscObject)subdm), &fvm);
115: PetscObjectSetName((PetscObject)fvm, "testField");
116: PetscFVSetNumComponents(fvm, 1);
117: DMGetCoordinateDim(subdm, &cdim);
118: PetscFVSetSpatialDimension(fvm, cdim);
119: PetscFVSetFromOptions(fvm);
121: DMAddField(subdm, NULL, (PetscObject)fvm);
122: PetscFVDestroy(&fvm);
123: DMCreateDS(subdm);
125: DMCreateGlobalVector(subdm, &gv);
126: PetscObjectSetName((PetscObject)gv, "potential");
127: VecSet(gv, 1.);
128: VecViewFromOptions(gv, NULL, "-vec_view");
129: VecDestroy(&gv);
130: DMDestroy(&subdm);
131: return 0;
132: }
134: int main(int argc, char **argv)
135: {
136: DM dm, subdm;
137: PetscReal height = -1.0;
138: PetscBool volume = PETSC_TRUE, domain = PETSC_FALSE;
141: PetscInitialize(&argc, &argv, NULL, help);
142: PetscOptionsGetReal(NULL, NULL, "-height", &height, NULL);
143: PetscOptionsGetBool(NULL, NULL, "-volume", &volume, NULL);
144: PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL);
146: CreateMesh(PETSC_COMM_WORLD, &dm);
147: if (volume) {
148: CreateVolumeSubmesh(dm, domain, PETSC_TRUE, height, &subdm);
149: DMSetFromOptions(subdm);
150: DMDestroy(&subdm);
151: CreateVolumeSubmesh(dm, domain, PETSC_FALSE, height, &subdm);
152: DMSetFromOptions(subdm);
153: DMDestroy(&subdm);
154: } else {
155: TestBoundaryField(dm);
156: }
157: DMDestroy(&dm);
158: PetscFinalize();
159: return 0;
160: }
162: /*TEST
164: test:
165: suffix: 0
166: requires: triangle
167: args: -dm_coord_space 0 -sub_dm_plex_check_all \
168: -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view
170: # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes
171: testset:
172: nsize: 3
173: requires: parmetis
174: args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \
175: -sub_dm_plex_check_all -dm_view -sub_dm_view
177: test:
178: suffix: 1
180: test:
181: suffix: 2
182: args: -domain
184: # This set tests that global numberings can be made when some strata are missing on a process
185: testset:
186: nsize: 3
187: requires: hdf5
188: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
189: -sub_dm_plex_check_all -sub_dm_view hdf5:subdm.h5
191: test:
192: suffix: 3
193: args: -domain -height 0.625
195: # This test checks whether filter can extract a lower-dimensional manifold and output a field on it
196: testset:
197: args: -volume 0 -dm_plex_simplex 0 -sub_dm_view hdf5:subdm.h5 -vec_view hdf5:subdm.h5::append
198: requires: hdf5
200: test:
201: suffix: surface_2d
202: args: -dm_plex_box_faces 5,5
204: test:
205: suffix: surface_3d
206: args: -dm_plex_box_faces 5,5,5
208: # This test checks that if cells are present in the SF, the dm is marked as having an overlap
209: test:
210: suffix: surface_2d_2
211: nsize: 3
212: args: -dm_plex_box_faces 5,5 -domain -height 0.625
214: TEST*/