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