xref: /petsc/src/dm/label/tutorials/ex1.c (revision df4cd43f92eaa320656440c40edb1046daee8f75)
1 static char help[] = "Tests DMLabel operations.\n\n";
2 
3 #include <petscdm.h>
4 #include <petscdmplex.h>
5 #include <petscdmplextransform.h>
6 
7 PetscErrorCode ViewLabels(DM dm, PetscViewer viewer)
8 {
9   DMLabel     label;
10   const char *labelName, *typeName;
11   PetscInt    numLabels, l;
12 
13   PetscFunctionBegin;
14   /* query the number and name of labels*/
15   PetscCall(DMGetNumLabels(dm, &numLabels));
16   PetscCall(PetscViewerASCIIPrintf(viewer, "Number of labels: %" PetscInt_FMT "\n", numLabels));
17   for (l = 0; l < numLabels; ++l) {
18     IS labelIS, tmpIS;
19 
20     PetscCall(DMGetLabelName(dm, l, &labelName));
21     PetscCall(DMGetLabel(dm, labelName, &label));
22     PetscCall(DMLabelGetType(label, &typeName));
23     PetscCall(PetscViewerASCIIPrintf(viewer, "Label %" PetscInt_FMT ": name: %s type: %s\n", l, labelName, typeName));
24     PetscCall(PetscViewerASCIIPrintf(viewer, "IS of values\n"));
25     PetscCall(DMLabelGetValueIS(label, &labelIS));
26     PetscCall(ISOnComm(labelIS, PetscObjectComm((PetscObject)viewer), PETSC_USE_POINTER, &tmpIS));
27     PetscCall(PetscViewerASCIIPushTab(viewer));
28     PetscCall(ISView(tmpIS, viewer));
29     PetscCall(PetscViewerASCIIPopTab(viewer));
30     PetscCall(ISDestroy(&tmpIS));
31     PetscCall(ISDestroy(&labelIS));
32     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
33   }
34   /* Making sure that string literals work */
35   PetscCall(PetscViewerASCIIPrintf(viewer, "\n\nCell Set label IS\n"));
36   PetscCall(DMGetLabel(dm, "Cell Sets", &label));
37   if (label) {
38     IS labelIS, tmpIS;
39 
40     PetscCall(DMLabelGetValueIS(label, &labelIS));
41     PetscCall(ISOnComm(labelIS, PetscObjectComm((PetscObject)viewer), PETSC_USE_POINTER, &tmpIS));
42     PetscCall(ISView(tmpIS, viewer));
43     PetscCall(ISDestroy(&tmpIS));
44     PetscCall(ISDestroy(&labelIS));
45   }
46   PetscFunctionReturn(PETSC_SUCCESS);
47 }
48 
49 PetscErrorCode CheckLabelsSame(DMLabel label0, DMLabel label1)
50 {
51   const char *name0, *name1;
52   PetscBool   same;
53   char       *msg;
54 
55   PetscFunctionBegin;
56   PetscCall(PetscObjectGetName((PetscObject)label0, &name0));
57   PetscCall(PetscObjectGetName((PetscObject)label1, &name1));
58   PetscCall(DMLabelCompare(PETSC_COMM_WORLD, label0, label1, &same, &msg));
59   PetscCheck(same == (PetscBool)!msg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "DMLabelCompare returns inconsistent same=%d msg=\"%s\"", same, msg);
60   PetscCheck(same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Labels \"%s\" and \"%s\" should not differ! Message:\n%s", name0, name1, msg);
61   /* Test passing NULL, must not fail */
62   PetscCall(DMLabelCompare(PETSC_COMM_WORLD, label0, label1, NULL, NULL));
63   PetscCall(PetscFree(msg));
64   PetscFunctionReturn(PETSC_SUCCESS);
65 }
66 
67 PetscErrorCode CheckLabelsNotSame(DMLabel label0, DMLabel label1)
68 {
69   const char *name0, *name1;
70   PetscBool   same;
71   char       *msg;
72 
73   PetscFunctionBegin;
74   PetscCall(PetscObjectGetName((PetscObject)label0, &name0));
75   PetscCall(PetscObjectGetName((PetscObject)label1, &name1));
76   PetscCall(DMLabelCompare(PETSC_COMM_WORLD, label0, label1, &same, &msg));
77   PetscCheck(same == (PetscBool)!msg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "DMLabelCompare returns inconsistent same=%d msg=\"%s\"", same, msg);
78   PetscCheck(!same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Labels \"%s\" and \"%s\" should differ!", name0, name1);
79   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Compare label \"%s\" with \"%s\": %s\n", name0, name1, msg));
80   PetscCall(PetscFree(msg));
81   PetscFunctionReturn(PETSC_SUCCESS);
82 }
83 
84 PetscErrorCode CheckDMLabelsSame(DM dm0, DM dm1)
85 {
86   const char *name0, *name1;
87   PetscBool   same;
88   char       *msg;
89 
90   PetscFunctionBegin;
91   PetscCall(PetscObjectGetName((PetscObject)dm0, &name0));
92   PetscCall(PetscObjectGetName((PetscObject)dm1, &name1));
93   PetscCall(DMCompareLabels(dm0, dm1, &same, &msg));
94   PetscCheck(same == (PetscBool)!msg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "DMCompareLabels returns inconsistent same=%d msg=\"%s\"", same, msg);
95   PetscCheck(same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Labels of DMs \"%s\" and \"%s\" should not differ! Message:\n%s", name0, name1, msg);
96   /* Test passing NULL, must not fail */
97   PetscCall(DMCompareLabels(dm0, dm1, NULL, NULL));
98   PetscCall(PetscFree(msg));
99   PetscFunctionReturn(PETSC_SUCCESS);
100 }
101 
102 PetscErrorCode CheckDMLabelsNotSame(DM dm0, DM dm1)
103 {
104   const char *name0, *name1;
105   PetscBool   same;
106   char       *msg;
107 
108   PetscFunctionBegin;
109   PetscCall(PetscObjectGetName((PetscObject)dm0, &name0));
110   PetscCall(PetscObjectGetName((PetscObject)dm1, &name1));
111   PetscCall(DMCompareLabels(dm0, dm1, &same, &msg));
112   PetscCheck(same == (PetscBool)!msg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "DMCompareLabels returns inconsistent same=%d msg=\"%s\"", same, msg);
113   PetscCheck(!same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Labels of DMs \"%s\" and \"%s\" should differ!", name0, name1);
114   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Labels of DMs \"%s\" and \"%s\" differ: %s\n", name0, name1, msg));
115   PetscCall(PetscFree(msg));
116   PetscFunctionReturn(PETSC_SUCCESS);
117 }
118 
119 PetscErrorCode CreateMesh(const char name[], DM *newdm)
120 {
121   DM        dm, dmDist;
122   char      filename[PETSC_MAX_PATH_LEN] = "";
123   PetscBool interpolate                  = PETSC_FALSE;
124 
125   PetscFunctionBegin;
126   /* initialize and get options */
127   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "DMLabel ex1 Options", "DMLabel");
128   PetscCall(PetscOptionsString("-i", "filename to read", "ex1.c", filename, filename, sizeof(filename), NULL));
129   PetscCall(PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex1.c", interpolate, &interpolate, NULL));
130   PetscOptionsEnd();
131 
132   /* create and distribute DM */
133   PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, "ex1_plex", interpolate, &dm));
134   PetscCall(DMPlexDistribute(dm, 0, NULL, &dmDist));
135   if (dmDist) {
136     PetscCall(DMDestroy(&dm));
137     dm = dmDist;
138   }
139   PetscCall(DMSetFromOptions(dm));
140   PetscCall(PetscObjectSetName((PetscObject)dm, name));
141   *newdm = dm;
142   PetscFunctionReturn(PETSC_SUCCESS);
143 }
144 
145 static PetscErrorCode TestEphemeralLabels(DM dm)
146 {
147   DMPlexTransform tr;
148   DM              tdm;
149   DMLabel         label, labelTmp;
150 
151   PetscFunctionBeginUser;
152   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
153   PetscCall(PetscObjectSetName((PetscObject)tr, "Transform"));
154   PetscCall(DMPlexTransformSetDM(tr, dm));
155   PetscCall(DMPlexTransformSetFromOptions(tr));
156   PetscCall(DMPlexTransformSetUp(tr));
157 
158   PetscCall(DMPlexCreateEphemeral(tr, &tdm));
159   PetscCall(DMPlexTransformDestroy(&tr));
160   PetscCall(PetscObjectSetName((PetscObject)tdm, "Ephemeral Mesh"));
161   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tdm, "eph_"));
162 
163   PetscCall(DMGetLabel(tdm, "marker", &label));
164   PetscCall(DMLabelDuplicate(label, &labelTmp));
165   PetscCall(CheckLabelsSame(label, labelTmp));
166   PetscCall(DMLabelDestroy(&labelTmp));
167   PetscCall(DMDestroy(&tdm));
168   PetscFunctionReturn(PETSC_SUCCESS);
169 }
170 
171 int main(int argc, char **argv)
172 {
173   DM dm;
174 
175   PetscFunctionBeginUser;
176   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
177   PetscCall(CreateMesh("plex0", &dm));
178   /* add custom labels to test adding/removal */
179   {
180     DMLabel  label0, label1, label2, label3;
181     PetscInt p, pStart, pEnd;
182     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
183     /* create label in DM and get from DM */
184     PetscCall(DMCreateLabel(dm, "label0"));
185     PetscCall(DMGetLabel(dm, "label0", &label0));
186     /* alternative: create standalone label and add to DM; needs to be destroyed */
187     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "label1", &label1));
188     PetscCall(DMAddLabel(dm, label1));
189 
190     pEnd = PetscMin(pEnd, pStart + 5);
191     for (p = pStart; p < pEnd; p++) {
192       PetscCall(DMLabelSetValue(label0, p, 1));
193       PetscCall(DMLabelSetValue(label1, p, 2));
194     }
195     /* duplicate label */
196     PetscCall(DMLabelDuplicate(label0, &label2));
197     PetscCall(DMLabelDuplicate(label1, &label3));
198     PetscCall(PetscObjectSetName((PetscObject)label2, "label2"));
199     PetscCall(PetscObjectSetName((PetscObject)label3, "label3"));
200     PetscCall(DMAddLabel(dm, label2));
201     PetscCall(DMAddLabel(dm, label3));
202     /* remove the labels in this scope */
203     PetscCall(DMLabelDestroy(&label1));
204     PetscCall(DMLabelDestroy(&label2));
205     PetscCall(DMLabelDestroy(&label3));
206   }
207 
208   PetscCall(ViewLabels(dm, PETSC_VIEWER_STDOUT_WORLD));
209 
210   /* do label perturbations and comparisons */
211   {
212     DMLabel  label0, label1, label2, label3;
213     PetscInt val;
214     PetscInt p, pStart, pEnd;
215 
216     PetscCall(DMGetLabel(dm, "label0", &label0));
217     PetscCall(DMGetLabel(dm, "label1", &label1));
218     PetscCall(DMGetLabel(dm, "label2", &label2));
219     PetscCall(DMGetLabel(dm, "label3", &label3));
220 
221     PetscCall(CheckLabelsNotSame(label0, label1));
222     PetscCall(CheckLabelsSame(label0, label2));
223     PetscCall(CheckLabelsSame(label1, label3));
224 
225     PetscCall(DMLabelGetDefaultValue(label1, &val));
226     PetscCall(DMLabelSetDefaultValue(label1, 333));
227     PetscCall(CheckLabelsNotSame(label1, label3));
228     PetscCall(DMLabelSetDefaultValue(label1, val));
229     PetscCall(CheckLabelsSame(label1, label3));
230 
231     PetscCall(DMLabelGetBounds(label1, &pStart, &pEnd));
232 
233     for (p = pStart; p < pEnd; p++) {
234       PetscCall(DMLabelGetValue(label1, p, &val));
235       // This is weird. Perhaps we should not need to call DMLabelClearValue()
236       PetscCall(DMLabelClearValue(label1, p, val));
237       val++;
238       PetscCall(DMLabelSetValue(label1, p, val));
239     }
240     PetscCall(CheckLabelsNotSame(label1, label3));
241     for (p = pStart; p < pEnd; p++) {
242       PetscCall(DMLabelGetValue(label1, p, &val));
243       // This is weird. Perhaps we should not need to call DMLabelClearValue()
244       PetscCall(DMLabelClearValue(label1, p, val));
245       val--;
246       PetscCall(DMLabelSetValue(label1, p, val));
247     }
248     PetscCall(CheckLabelsSame(label1, label3));
249 
250     PetscCall(DMLabelGetValue(label3, pEnd - 1, &val));
251     PetscCall(DMLabelSetValue(label3, pEnd, val));
252     PetscCall(CheckLabelsNotSame(label1, label3));
253     // This is weird. Perhaps we should not need to call DMLabelClearValue()
254     PetscCall(DMLabelClearValue(label3, pEnd, val));
255     PetscCall(CheckLabelsSame(label1, label3));
256   }
257 
258   {
259     DM       dm1;
260     DMLabel  label02, label12;
261     PetscInt p = 0, val;
262 
263     PetscCall(CreateMesh("plex1", &dm1));
264     PetscCall(CheckDMLabelsNotSame(dm, dm1));
265 
266     PetscCall(DMCopyLabels(dm, dm1, PETSC_OWN_POINTER, PETSC_FALSE, DM_COPY_LABELS_REPLACE));
267     PetscCall(CheckDMLabelsSame(dm, dm1));
268 
269     PetscCall(DMCopyLabels(dm, dm1, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_REPLACE));
270     PetscCall(DMGetLabel(dm, "label2", &label02));
271     PetscCall(DMGetLabel(dm1, "label2", &label12));
272     PetscCall(CheckLabelsSame(label02, label12));
273 
274     PetscCall(DMLabelGetValue(label12, p, &val));
275     // This is weird. Perhaps we should not need to call DMLabelClearValue()
276     PetscCall(DMLabelClearValue(label12, p, val));
277     PetscCall(DMLabelSetValue(label12, p, val + 1));
278     PetscCall(CheckLabelsNotSame(label02, label12));
279     PetscCall(CheckDMLabelsNotSame(dm, dm1));
280 
281     // This is weird. Perhaps we should not need to call DMLabelClearValue()
282     PetscCall(DMLabelClearValue(label12, p, val + 1));
283     PetscCall(DMLabelSetValue(label12, p, val));
284     PetscCall(CheckLabelsSame(label02, label12));
285     PetscCall(CheckDMLabelsSame(dm, dm1));
286 
287     PetscCall(PetscObjectSetName((PetscObject)label12, "label12"));
288     PetscCall(CheckDMLabelsNotSame(dm, dm1));
289     PetscCall(PetscObjectSetName((PetscObject)label12, "label2"));
290     PetscCall(CheckDMLabelsSame(dm, dm1));
291 
292     PetscCall(DMDestroy(&dm1));
293   }
294   // Test adding strata and filtering
295   {
296     DMLabel  labelA, labelB;
297     IS       valueIS;
298     PetscInt pStart, pEnd, lStart = 5, lEnd = 19;
299 
300     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
301     PetscCall(DMCreateLabel(dm, "labelA"));
302     PetscCall(DMCreateLabel(dm, "labelB"));
303     PetscCall(DMGetLabel(dm, "labelA", &labelA));
304     PetscCall(DMGetLabel(dm, "labelB", &labelB));
305     for (PetscInt p = pStart; p < pEnd; ++p) {
306       if (p < lStart || p >= lEnd) continue;
307       if (p % 2) PetscCall(DMLabelSetValue(labelA, p, 19));
308       else PetscCall(DMLabelSetValue(labelA, p, 17));
309     }
310     PetscCall(DMLabelGetValueIS(labelA, &valueIS));
311     PetscCall(DMLabelAddStrataIS(labelA, valueIS));
312     PetscCall(ISDestroy(&valueIS));
313     for (PetscInt p = pStart; p < pEnd; ++p) {
314       if (p % 2) PetscCall(DMLabelSetValue(labelB, p, 19));
315       else PetscCall(DMLabelSetValue(labelB, p, 17));
316     }
317     PetscCall(DMLabelFilter(labelB, lStart, lEnd));
318     PetscCall(CheckLabelsSame(labelA, labelB));
319     PetscCall(DMRemoveLabel(dm, "labelA", NULL));
320     PetscCall(DMRemoveLabel(dm, "labelB", NULL));
321   }
322 
323   /* remove label0 and label1 just to test manual removal; let label3 be removed automatically by DMDestroy() */
324   {
325     DMLabel label0, label1, label2;
326     PetscCall(DMGetLabel(dm, "label0", &label0));
327     PetscCall(DMGetLabel(dm, "label1", &label1));
328     PetscCheck(label0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label0 must not be NULL now");
329     PetscCheck(label1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label1 must not be NULL now");
330     PetscCall(DMRemoveLabel(dm, "label1", NULL));
331     PetscCall(DMRemoveLabel(dm, "label2", &label2));
332     PetscCall(DMRemoveLabelBySelf(dm, &label0, PETSC_TRUE));
333     PetscCall(DMGetLabel(dm, "label0", &label0));
334     PetscCall(DMGetLabel(dm, "label1", &label1));
335     PetscCheck(!label0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label0 must be NULL now");
336     PetscCheck(!label1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label1 must be NULL now");
337     PetscCheck(label2, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label2 must not be NULL now");
338     PetscCall(DMRemoveLabelBySelf(dm, &label2, PETSC_FALSE)); /* this should do nothing */
339     PetscCheck(label2, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label2 must not be NULL now");
340     PetscCall(DMLabelDestroy(&label2));
341     PetscCall(DMGetLabel(dm, "label2", &label2));
342     PetscCheck(!label2, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "label2 must be NULL now");
343   }
344 
345   PetscCall(TestEphemeralLabels(dm));
346 
347   PetscCall(DMDestroy(&dm));
348   PetscCall(PetscFinalize());
349   return 0;
350 }
351 
352 /*TEST
353 
354   test:
355     suffix: 0
356     nsize: {{1 2}separate output}
357     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -interpolate
358     requires: exodusii
359 
360 TEST*/
361