xref: /petsc/src/dm/impls/plex/tests/ex11.c (revision d71ae5a4db6382e7f06317b8d368875286fe9008)
1 static char help[] = "Tests for DMLabel\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petsc/private/dmimpl.h>
5 
6 static PetscErrorCode TestInsertion()
7 {
8   DMLabel        label, label2;
9   const PetscInt values[5] = {0, 3, 4, -1, 176}, N = 10000;
10   PetscInt       i, v;
11 
12   PetscFunctionBegin;
13   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Test Label", &label));
14   PetscCall(DMLabelSetDefaultValue(label, -100));
15   for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, i, values[i % 5]));
16   /* Test get in hash mode */
17   for (i = 0; i < N; ++i) {
18     PetscInt val;
19 
20     PetscCall(DMLabelGetValue(label, i, &val));
21     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " for point %" PetscInt_FMT " should be %" PetscInt_FMT, val, i, values[i % 5]);
22   }
23   /* Test stratum */
24   for (v = 0; v < 5; ++v) {
25     IS              stratum;
26     const PetscInt *points;
27     PetscInt        n;
28 
29     PetscCall(DMLabelGetStratumIS(label, values[v], &stratum));
30     PetscCheck(stratum, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Stratum %" PetscInt_FMT " is empty!", v);
31     PetscCall(ISGetIndices(stratum, &points));
32     PetscCall(ISGetLocalSize(stratum, &n));
33     for (i = 0; i < n; ++i) PetscCheck(points[i] == i * 5 + v, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " should be %" PetscInt_FMT, points[i], i * 5 + v);
34     PetscCall(ISRestoreIndices(stratum, &points));
35     PetscCall(ISDestroy(&stratum));
36   }
37   /* Test get in array mode */
38   for (i = 0; i < N; ++i) {
39     PetscInt val;
40 
41     PetscCall(DMLabelGetValue(label, i, &val));
42     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]);
43   }
44   /* Test Duplicate */
45   PetscCall(DMLabelDuplicate(label, &label2));
46   for (i = 0; i < N; ++i) {
47     PetscInt val;
48 
49     PetscCall(DMLabelGetValue(label2, i, &val));
50     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]);
51   }
52   PetscCall(DMLabelDestroy(&label2));
53   PetscCall(DMLabelDestroy(&label));
54   PetscFunctionReturn(0);
55 }
56 
57 static PetscErrorCode TestEmptyStrata(MPI_Comm comm)
58 {
59   DM               dm, dmDist;
60   PetscPartitioner part;
61   PetscInt         c0[6]  = {2, 3, 6, 7, 9, 11};
62   PetscInt         c1[6]  = {4, 5, 7, 8, 10, 12};
63   PetscInt         c2[4]  = {13, 15, 19, 21};
64   PetscInt         c3[4]  = {14, 16, 20, 22};
65   PetscInt         c4[4]  = {15, 17, 21, 23};
66   PetscInt         c5[4]  = {16, 18, 22, 24};
67   PetscInt         c6[4]  = {13, 14, 19, 20};
68   PetscInt         c7[4]  = {15, 16, 21, 22};
69   PetscInt         c8[4]  = {17, 18, 23, 24};
70   PetscInt         c9[4]  = {13, 14, 15, 16};
71   PetscInt         c10[4] = {15, 16, 17, 18};
72   PetscInt         c11[4] = {19, 20, 21, 22};
73   PetscInt         c12[4] = {21, 22, 23, 24};
74   PetscInt         dim    = 3;
75   PetscMPIInt      rank;
76 
77   PetscFunctionBegin;
78   PetscCallMPI(MPI_Comm_rank(comm, &rank));
79   /* A 3D box with two adjacent cells, sharing one face and four vertices */
80   PetscCall(DMCreate(comm, &dm));
81   PetscCall(DMSetType(dm, DMPLEX));
82   PetscCall(DMSetDimension(dm, dim));
83   if (rank == 0) {
84     PetscCall(DMPlexSetChart(dm, 0, 25));
85     PetscCall(DMPlexSetConeSize(dm, 0, 6));
86     PetscCall(DMPlexSetConeSize(dm, 1, 6));
87     PetscCall(DMPlexSetConeSize(dm, 2, 4));
88     PetscCall(DMPlexSetConeSize(dm, 3, 4));
89     PetscCall(DMPlexSetConeSize(dm, 4, 4));
90     PetscCall(DMPlexSetConeSize(dm, 5, 4));
91     PetscCall(DMPlexSetConeSize(dm, 6, 4));
92     PetscCall(DMPlexSetConeSize(dm, 7, 4));
93     PetscCall(DMPlexSetConeSize(dm, 8, 4));
94     PetscCall(DMPlexSetConeSize(dm, 9, 4));
95     PetscCall(DMPlexSetConeSize(dm, 10, 4));
96     PetscCall(DMPlexSetConeSize(dm, 11, 4));
97     PetscCall(DMPlexSetConeSize(dm, 12, 4));
98   }
99   PetscCall(DMSetUp(dm));
100   if (rank == 0) {
101     PetscCall(DMPlexSetCone(dm, 0, c0));
102     PetscCall(DMPlexSetCone(dm, 1, c1));
103     PetscCall(DMPlexSetCone(dm, 2, c2));
104     PetscCall(DMPlexSetCone(dm, 3, c3));
105     PetscCall(DMPlexSetCone(dm, 4, c4));
106     PetscCall(DMPlexSetCone(dm, 5, c5));
107     PetscCall(DMPlexSetCone(dm, 6, c6));
108     PetscCall(DMPlexSetCone(dm, 7, c7));
109     PetscCall(DMPlexSetCone(dm, 8, c8));
110     PetscCall(DMPlexSetCone(dm, 9, c9));
111     PetscCall(DMPlexSetCone(dm, 10, c10));
112     PetscCall(DMPlexSetCone(dm, 11, c11));
113     PetscCall(DMPlexSetCone(dm, 12, c12));
114   }
115   PetscCall(DMPlexSymmetrize(dm));
116   /* Create a user managed depth label, so that we can leave out edges */
117   {
118     DMLabel  label;
119     PetscInt numValues, maxValues = 0, v;
120 
121     PetscCall(DMCreateLabel(dm, "depth"));
122     PetscCall(DMPlexGetDepthLabel(dm, &label));
123     if (rank == 0) {
124       PetscInt i;
125 
126       for (i = 0; i < 25; ++i) {
127         if (i < 2) PetscCall(DMLabelSetValue(label, i, 3));
128         else if (i < 13) PetscCall(DMLabelSetValue(label, i, 2));
129         else {
130           if (i == 13) PetscCall(DMLabelAddStratum(label, 1));
131           PetscCall(DMLabelSetValue(label, i, 0));
132         }
133       }
134     }
135     PetscCall(DMLabelGetNumValues(label, &numValues));
136     PetscCallMPI(MPI_Allreduce(&numValues, &maxValues, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
137     for (v = numValues; v < maxValues; ++v) PetscCall(DMLabelAddStratum(label, v));
138   }
139   {
140     DMLabel label;
141     PetscCall(DMPlexGetDepthLabel(dm, &label));
142     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm)));
143   }
144   PetscCall(DMPlexGetPartitioner(dm, &part));
145   PetscCall(PetscPartitionerSetFromOptions(part));
146   PetscCall(DMPlexDistribute(dm, 1, NULL, &dmDist));
147   if (dmDist) {
148     PetscCall(DMDestroy(&dm));
149     dm = dmDist;
150   }
151   {
152     DMLabel label;
153     PetscCall(DMPlexGetDepthLabel(dm, &label));
154     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm)));
155   }
156   /* Create a cell vector */
157   {
158     Vec          v;
159     PetscSection s;
160     PetscInt     numComp[] = {1};
161     PetscInt     dof[]     = {0, 0, 0, 1};
162     PetscInt     N;
163 
164     PetscCall(DMSetNumFields(dm, 1));
165     PetscCall(DMPlexCreateSection(dm, NULL, numComp, dof, 0, NULL, NULL, NULL, NULL, &s));
166     PetscCall(DMSetLocalSection(dm, s));
167     PetscCall(PetscSectionDestroy(&s));
168     PetscCall(DMCreateGlobalVector(dm, &v));
169     PetscCall(VecGetSize(v, &N));
170     if (N != 2) {
171       PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_(comm)));
172       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "FAIL: Vector size %" PetscInt_FMT " != 2", N);
173     }
174     PetscCall(VecDestroy(&v));
175   }
176   PetscCall(DMDestroy(&dm));
177   PetscFunctionReturn(0);
178 }
179 
180 static PetscErrorCode TestDistribution(MPI_Comm comm)
181 {
182   DM               dm, dmDist;
183   PetscPartitioner part;
184   DMLabel          label;
185   char             filename[PETSC_MAX_PATH_LEN];
186   const char      *name    = "test label";
187   PetscInt         overlap = 0, cStart, cEnd, c;
188   PetscMPIInt      rank;
189   PetscBool        flg;
190 
191   PetscFunctionBegin;
192   PetscCallMPI(MPI_Comm_rank(comm, &rank));
193   PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg));
194   if (!flg) PetscFunctionReturn(0);
195   PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL));
196   PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm));
197   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
198   PetscCall(DMCreateLabel(dm, name));
199   PetscCall(DMGetLabel(dm, name, &label));
200   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
201   for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(label, c, c));
202   PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
203   PetscCall(DMPlexGetPartitioner(dm, &part));
204   PetscCall(PetscPartitionerSetFromOptions(part));
205   PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist));
206   if (dmDist) {
207     PetscCall(DMDestroy(&dm));
208     dm = dmDist;
209   }
210   PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
211   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
212   PetscCall(DMGetLabel(dm, name, &label));
213   PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
214   PetscCall(DMDestroy(&dm));
215   PetscFunctionReturn(0);
216 }
217 
218 static PetscErrorCode TestUniversalLabel(MPI_Comm comm)
219 {
220   DM               dm1, dm2;
221   DMLabel          bd1, bd2, ulabel;
222   DMUniversalLabel universal;
223   PetscInt         pStart, pEnd, p;
224   PetscBool        run = PETSC_FALSE, notFile;
225 
226   PetscFunctionBeginUser;
227   PetscCall(PetscOptionsGetBool(NULL, NULL, "-universal", &run, NULL));
228   if (!run) PetscFunctionReturn(0);
229 
230   char      filename[PETSC_MAX_PATH_LEN];
231   PetscBool flg;
232 
233   PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg));
234   if (flg) {
235     PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm1));
236   } else {
237     PetscCall(DMCreate(comm, &dm1));
238     PetscCall(DMSetType(dm1, DMPLEX));
239     PetscCall(DMSetFromOptions(dm1));
240   }
241   PetscCall(DMHasLabel(dm1, "marker", &notFile));
242   if (notFile) {
243     PetscCall(DMCreateLabel(dm1, "Boundary Faces"));
244     PetscCall(DMGetLabel(dm1, "Boundary Faces", &bd1));
245     PetscCall(DMPlexMarkBoundaryFaces(dm1, 13, bd1));
246     PetscCall(DMCreateLabel(dm1, "Boundary"));
247     PetscCall(DMGetLabel(dm1, "Boundary", &bd2));
248     PetscCall(DMPlexMarkBoundaryFaces(dm1, 121, bd2));
249     PetscCall(DMPlexLabelComplete(dm1, bd2));
250   }
251   PetscCall(PetscObjectSetName((PetscObject)dm1, "First Mesh"));
252   PetscCall(DMViewFromOptions(dm1, NULL, "-dm_view"));
253 
254   PetscCall(DMUniversalLabelCreate(dm1, &universal));
255   PetscCall(DMUniversalLabelGetLabel(universal, &ulabel));
256   PetscCall(PetscObjectViewFromOptions((PetscObject)ulabel, NULL, "-universal_view"));
257 
258   if (!notFile) {
259     PetscInt Nl, l;
260 
261     PetscCall(DMClone(dm1, &dm2));
262     PetscCall(DMGetNumLabels(dm2, &Nl));
263     for (l = Nl - 1; l >= 0; --l) {
264       PetscBool   isdepth, iscelltype;
265       const char *name;
266 
267       PetscCall(DMGetLabelName(dm2, l, &name));
268       PetscCall(PetscStrncmp(name, "depth", 6, &isdepth));
269       PetscCall(PetscStrncmp(name, "celltype", 9, &iscelltype));
270       if (!isdepth && !iscelltype) PetscCall(DMRemoveLabel(dm2, name, NULL));
271     }
272   } else {
273     PetscCall(DMCreate(comm, &dm2));
274     PetscCall(DMSetType(dm2, DMPLEX));
275     PetscCall(DMSetFromOptions(dm2));
276   }
277   PetscCall(PetscObjectSetName((PetscObject)dm2, "Second Mesh"));
278   PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, dm2));
279   PetscCall(DMPlexGetChart(dm2, &pStart, &pEnd));
280   for (p = pStart; p < pEnd; ++p) {
281     PetscInt val;
282 
283     PetscCall(DMLabelGetValue(ulabel, p, &val));
284     if (val < 0) continue;
285     PetscCall(DMUniversalLabelSetLabelValue(universal, dm2, PETSC_TRUE, p, val));
286   }
287   PetscCall(DMViewFromOptions(dm2, NULL, "-dm_view"));
288 
289   PetscCall(DMUniversalLabelDestroy(&universal));
290   PetscCall(DMDestroy(&dm1));
291   PetscCall(DMDestroy(&dm2));
292   PetscFunctionReturn(0);
293 }
294 
295 int main(int argc, char **argv)
296 {
297   PetscFunctionBeginUser;
298   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
299   /*PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));*/
300   PetscCall(TestInsertion());
301   PetscCall(TestEmptyStrata(PETSC_COMM_WORLD));
302   PetscCall(TestDistribution(PETSC_COMM_WORLD));
303   PetscCall(TestUniversalLabel(PETSC_COMM_WORLD));
304   PetscCall(PetscFinalize());
305   return 0;
306 }
307 
308 /*TEST
309 
310   test:
311     suffix: 0
312     requires: triangle
313   test:
314     suffix: 1
315     requires: triangle
316     nsize: 2
317     args: -petscpartitioner_type simple
318 
319   testset:
320     suffix: gmsh
321     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square.msh -petscpartitioner_type simple
322     test:
323       suffix: 1
324       nsize: 1
325     test:
326       suffix: 2
327       nsize: 2
328 
329   testset:
330     suffix: exodusii
331     requires: exodusii
332     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/2Dgrd.exo -petscpartitioner_type simple
333     test:
334       suffix: 1
335       nsize: 1
336     test:
337       suffix: 2
338       nsize: 2
339 
340   test:
341     suffix: univ
342     requires: triangle
343     args: -universal -dm_view -universal_view
344 
345   test:
346     # Note that the labels differ because we have multiply-marked some points during EGADS creation
347     suffix: univ_egads_sphere
348     requires: egads
349     args: -universal -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_view -universal_view
350 
351   test:
352     # Note that the labels differ because we have multiply-marked some points during EGADS creation
353     suffix: univ_egads_ball
354     requires: egads ctetgen
355     args: -universal -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_view -universal_view
356 
357 TEST*/
358