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