1 static char help[] = "Tests for creation of submeshes\n\n";
2
3 #include <petscdmplex.h>
4
CreateMesh(MPI_Comm comm,DM * dm)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 }
14
15 // Label half of the cells
CreateHalfCellsLabel(DM dm,PetscBool lower,DMLabel * label)16 static PetscErrorCode CreateHalfCellsLabel(DM dm, PetscBool lower, DMLabel *label)
17 {
18 PetscInt cStart, cEnd, cStartSub, cEndSub;
19
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 }
36
37 // Label everything on the right half of the domain
CreateHalfDomainLabel(DM dm,PetscBool lower,PetscReal height,DMLabel * label)38 static PetscErrorCode CreateHalfDomainLabel(DM dm, PetscBool lower, PetscReal height, DMLabel *label)
39 {
40 PetscReal centroid[3];
41 PetscInt cStart, cEnd, cdim;
42
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 }
62
63 // Create a line of faces at a given x value
CreateLineLabel(DM dm,PetscReal x,DMLabel * label)64 static PetscErrorCode CreateLineLabel(DM dm, PetscReal x, DMLabel *label)
65 {
66 PetscReal centroid[3];
67 PetscInt fStart, fEnd;
68
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 }
81
CreateVolumeSubmesh(DM dm,PetscBool domain,PetscBool lower,PetscReal height,DM * subdm)82 static PetscErrorCode CreateVolumeSubmesh(DM dm, PetscBool domain, PetscBool lower, PetscReal height, DM *subdm)
83 {
84 DMLabel label, map;
85
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, PetscObjectComm((PetscObject)dm), 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 }
97
TestBoundaryField(DM dm)98 static PetscErrorCode TestBoundaryField(DM dm)
99 {
100 DM subdm;
101 DMLabel label, map;
102 PetscFV fvm;
103 Vec gv;
104 PetscInt cdim;
105
106 PetscFunctionBeginUser;
107 PetscCall(CreateLineLabel(dm, 0.5, &label));
108 PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)dm), 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"));
114
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));
121
122 PetscCall(DMAddField(subdm, NULL, (PetscObject)fvm));
123 PetscCall(PetscFVDestroy(&fvm));
124 PetscCall(DMCreateDS(subdm));
125
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 }
134
main(int argc,char ** argv)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;
140
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));
146
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 }
162
163 /*TEST
164
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
170
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
176
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
183
184 test:
185 suffix: 1
186
187 test:
188 suffix: 2
189 args: -domain
190
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
197
198 test:
199 suffix: 3
200 args: -domain -height 0.625
201 output_file: output/empty.out
202
203 # This set tests that global numberings can be made when some strata are missing on a process
204 testset:
205 nsize: 3
206 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
207 -sub_dm_plex_check_all -sub_dm_view {{vtk:subdm.vtu :subdm.txt :subdm_d.txt:ascii_info_detail}}
208
209 test:
210 suffix: 3_vtk
211 args: -domain -height 0.625
212 output_file: output/empty.out
213
214 # This test checks whether filter can extract a lower-dimensional manifold and output a field on it
215 testset:
216 requires: hdf5
217 args: -volume 0 -dm_plex_simplex 0 -sub_dm_view hdf5:subdm.h5 -vec_view hdf5:subdm.h5::append
218 output_file: output/empty.out
219
220 test:
221 suffix: surface_2d
222 args: -dm_plex_box_faces 5,5
223
224 test:
225 suffix: surface_3d
226 args: -dm_plex_box_faces 5,5,5
227
228 # This test checks that if cells are present in the SF, the dm is marked as having an overlap
229 test:
230 suffix: surface_2d_2
231 nsize: 3
232 args: -dm_plex_box_faces 5,5 -domain -height 0.625
233
234 TEST*/
235