1 #include <petscdm.h>
2 #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
3 #include <petsc/private/sectionimpl.h> /*I "petscsection.h" I*/
4 #include <petscsf.h>
5 #include <petscsection.h>
6
7 /*@C
8 DMLabelCreate - Create a DMLabel object, which is a multimap
9
10 Collective
11
12 Input parameters:
13 + comm - The communicator, usually PETSC_COMM_SELF
14 - name - The label name
15
16 Output parameter:
17 . label - The DMLabel
18
19 Level: beginner
20
21 Notes:
22 The label name is actually usual PetscObject name.
23 One can get/set it with PetscObjectGetName()/PetscObjectSetName().
24
25 .seealso: DMLabelDestroy()
26 @*/
DMLabelCreate(MPI_Comm comm,const char name[],DMLabel * label)27 PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
28 {
29 PetscErrorCode ierr;
30
31 PetscFunctionBegin;
32 PetscValidPointer(label,2);
33 ierr = DMInitializePackage();CHKERRQ(ierr);
34
35 ierr = PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);CHKERRQ(ierr);
36
37 (*label)->numStrata = 0;
38 (*label)->defaultValue = -1;
39 (*label)->stratumValues = NULL;
40 (*label)->validIS = NULL;
41 (*label)->stratumSizes = NULL;
42 (*label)->points = NULL;
43 (*label)->ht = NULL;
44 (*label)->pStart = -1;
45 (*label)->pEnd = -1;
46 (*label)->bt = NULL;
47 ierr = PetscHMapICreate(&(*label)->hmap);CHKERRQ(ierr);
48 ierr = PetscObjectSetName((PetscObject) *label, name);CHKERRQ(ierr);
49 PetscFunctionReturn(0);
50 }
51
52 /*
53 DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
54
55 Not collective
56
57 Input parameter:
58 + label - The DMLabel
59 - v - The stratum value
60
61 Output parameter:
62 . label - The DMLabel with stratum in sorted list format
63
64 Level: developer
65
66 .seealso: DMLabelCreate()
67 */
DMLabelMakeValid_Private(DMLabel label,PetscInt v)68 static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
69 {
70 IS is;
71 PetscInt off = 0, *pointArray, p;
72 PetscErrorCode ierr;
73
74 if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
75 PetscFunctionBegin;
76 if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
77 ierr = PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);CHKERRQ(ierr);
78 ierr = PetscMalloc1(label->stratumSizes[v], &pointArray);CHKERRQ(ierr);
79 ierr = PetscHSetIGetElems(label->ht[v], &off, pointArray);CHKERRQ(ierr);
80 ierr = PetscHSetIClear(label->ht[v]);CHKERRQ(ierr);
81 ierr = PetscSortInt(label->stratumSizes[v], pointArray);CHKERRQ(ierr);
82 if (label->bt) {
83 for (p = 0; p < label->stratumSizes[v]; ++p) {
84 const PetscInt point = pointArray[p];
85 if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
86 ierr = PetscBTSet(label->bt, point - label->pStart);CHKERRQ(ierr);
87 }
88 }
89 if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
90 ierr = ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);CHKERRQ(ierr);
91 ierr = PetscFree(pointArray);CHKERRQ(ierr);
92 } else {
93 ierr = ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);CHKERRQ(ierr);
94 }
95 ierr = PetscObjectSetName((PetscObject) is, "indices");CHKERRQ(ierr);
96 label->points[v] = is;
97 label->validIS[v] = PETSC_TRUE;
98 ierr = PetscObjectStateIncrease((PetscObject) label);CHKERRQ(ierr);
99 PetscFunctionReturn(0);
100 }
101
102 /*
103 DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
104
105 Not collective
106
107 Input parameter:
108 . label - The DMLabel
109
110 Output parameter:
111 . label - The DMLabel with all strata in sorted list format
112
113 Level: developer
114
115 .seealso: DMLabelCreate()
116 */
DMLabelMakeAllValid_Private(DMLabel label)117 static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
118 {
119 PetscInt v;
120 PetscErrorCode ierr;
121
122 PetscFunctionBegin;
123 for (v = 0; v < label->numStrata; v++){
124 ierr = DMLabelMakeValid_Private(label, v);CHKERRQ(ierr);
125 }
126 PetscFunctionReturn(0);
127 }
128
129 /*
130 DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
131
132 Not collective
133
134 Input parameter:
135 + label - The DMLabel
136 - v - The stratum value
137
138 Output parameter:
139 . label - The DMLabel with stratum in hash format
140
141 Level: developer
142
143 .seealso: DMLabelCreate()
144 */
DMLabelMakeInvalid_Private(DMLabel label,PetscInt v)145 static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
146 {
147 PetscInt p;
148 const PetscInt *points;
149 PetscErrorCode ierr;
150
151 if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
152 PetscFunctionBegin;
153 if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
154 if (label->points[v]) {
155 ierr = ISGetIndices(label->points[v], &points);CHKERRQ(ierr);
156 for (p = 0; p < label->stratumSizes[v]; ++p) {
157 ierr = PetscHSetIAdd(label->ht[v], points[p]);CHKERRQ(ierr);
158 }
159 ierr = ISRestoreIndices(label->points[v],&points);CHKERRQ(ierr);
160 ierr = ISDestroy(&(label->points[v]));CHKERRQ(ierr);
161 }
162 label->validIS[v] = PETSC_FALSE;
163 PetscFunctionReturn(0);
164 }
165
166 #if !defined(DMLABEL_LOOKUP_THRESHOLD)
167 #define DMLABEL_LOOKUP_THRESHOLD 16
168 #endif
169
DMLabelLookupStratum(DMLabel label,PetscInt value,PetscInt * index)170 PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
171 {
172 PetscInt v;
173 PetscErrorCode ierr;
174
175 PetscFunctionBegin;
176 *index = -1;
177 if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
178 for (v = 0; v < label->numStrata; ++v)
179 if (label->stratumValues[v] == value) {*index = v; break;}
180 } else {
181 ierr = PetscHMapIGet(label->hmap, value, index);CHKERRQ(ierr);
182 }
183 if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
184 PetscInt len, loc = -1;
185 ierr = PetscHMapIGetSize(label->hmap, &len);CHKERRQ(ierr);
186 if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
187 if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
188 ierr = PetscHMapIGet(label->hmap, value, &loc);CHKERRQ(ierr);
189 } else {
190 for (v = 0; v < label->numStrata; ++v)
191 if (label->stratumValues[v] == value) {loc = v; break;}
192 }
193 if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
194 }
195 PetscFunctionReturn(0);
196 }
197
DMLabelNewStratum(DMLabel label,PetscInt value,PetscInt * index)198 PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
199 {
200 PetscInt v;
201 PetscInt *tmpV;
202 PetscInt *tmpS;
203 PetscHSetI *tmpH, ht;
204 IS *tmpP, is;
205 PetscBool *tmpB;
206 PetscHMapI hmap = label->hmap;
207 PetscErrorCode ierr;
208
209 PetscFunctionBegin;
210 v = label->numStrata;
211 tmpV = label->stratumValues;
212 tmpS = label->stratumSizes;
213 tmpH = label->ht;
214 tmpP = label->points;
215 tmpB = label->validIS;
216 { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
217 PetscInt *oldV = tmpV;
218 PetscInt *oldS = tmpS;
219 PetscHSetI *oldH = tmpH;
220 IS *oldP = tmpP;
221 PetscBool *oldB = tmpB;
222 ierr = PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);CHKERRQ(ierr);
223 ierr = PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);CHKERRQ(ierr);
224 ierr = PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);CHKERRQ(ierr);
225 ierr = PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);CHKERRQ(ierr);
226 ierr = PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);CHKERRQ(ierr);
227 ierr = PetscArraycpy(tmpV, oldV, v);CHKERRQ(ierr);
228 ierr = PetscArraycpy(tmpS, oldS, v);CHKERRQ(ierr);
229 ierr = PetscArraycpy(tmpH, oldH, v);CHKERRQ(ierr);
230 ierr = PetscArraycpy(tmpP, oldP, v);CHKERRQ(ierr);
231 ierr = PetscArraycpy(tmpB, oldB, v);CHKERRQ(ierr);
232 ierr = PetscFree(oldV);CHKERRQ(ierr);
233 ierr = PetscFree(oldS);CHKERRQ(ierr);
234 ierr = PetscFree(oldH);CHKERRQ(ierr);
235 ierr = PetscFree(oldP);CHKERRQ(ierr);
236 ierr = PetscFree(oldB);CHKERRQ(ierr);
237 }
238 label->numStrata = v+1;
239 label->stratumValues = tmpV;
240 label->stratumSizes = tmpS;
241 label->ht = tmpH;
242 label->points = tmpP;
243 label->validIS = tmpB;
244 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr);
245 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);CHKERRQ(ierr);
246 ierr = PetscHMapISet(hmap, value, v);CHKERRQ(ierr);
247 tmpV[v] = value;
248 tmpS[v] = 0;
249 tmpH[v] = ht;
250 tmpP[v] = is;
251 tmpB[v] = PETSC_TRUE;
252 ierr = PetscObjectStateIncrease((PetscObject) label);CHKERRQ(ierr);
253 *index = v;
254 PetscFunctionReturn(0);
255 }
256
DMLabelLookupAddStratum(DMLabel label,PetscInt value,PetscInt * index)257 PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
258 {
259 PetscErrorCode ierr;
260 PetscFunctionBegin;
261 ierr = DMLabelLookupStratum(label, value, index);CHKERRQ(ierr);
262 if (*index < 0) {ierr = DMLabelNewStratum(label, value, index);CHKERRQ(ierr);}
263 PetscFunctionReturn(0);
264 }
265
266 /*@
267 DMLabelAddStratum - Adds a new stratum value in a DMLabel
268
269 Input Parameter:
270 + label - The DMLabel
271 - value - The stratum value
272
273 Level: beginner
274
275 .seealso: DMLabelCreate(), DMLabelDestroy()
276 @*/
DMLabelAddStratum(DMLabel label,PetscInt value)277 PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
278 {
279 PetscInt v;
280 PetscErrorCode ierr;
281
282 PetscFunctionBegin;
283 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
284 ierr = DMLabelLookupAddStratum(label, value, &v);CHKERRQ(ierr);
285 PetscFunctionReturn(0);
286 }
287
288 /*@
289 DMLabelAddStrata - Adds new stratum values in a DMLabel
290
291 Not collective
292
293 Input Parameter:
294 + label - The DMLabel
295 . numStrata - The number of stratum values
296 - stratumValues - The stratum values
297
298 Level: beginner
299
300 .seealso: DMLabelCreate(), DMLabelDestroy()
301 @*/
DMLabelAddStrata(DMLabel label,PetscInt numStrata,const PetscInt stratumValues[])302 PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
303 {
304 PetscInt *values, v;
305 PetscErrorCode ierr;
306
307 PetscFunctionBegin;
308 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
309 if (numStrata) PetscValidIntPointer(stratumValues, 3);
310 ierr = PetscMalloc1(numStrata, &values);CHKERRQ(ierr);
311 ierr = PetscArraycpy(values, stratumValues, numStrata);CHKERRQ(ierr);
312 ierr = PetscSortRemoveDupsInt(&numStrata, values);CHKERRQ(ierr);
313 if (!label->numStrata) { /* Fast preallocation */
314 PetscInt *tmpV;
315 PetscInt *tmpS;
316 PetscHSetI *tmpH, ht;
317 IS *tmpP, is;
318 PetscBool *tmpB;
319 PetscHMapI hmap = label->hmap;
320
321 ierr = PetscMalloc1(numStrata, &tmpV);CHKERRQ(ierr);
322 ierr = PetscMalloc1(numStrata, &tmpS);CHKERRQ(ierr);
323 ierr = PetscMalloc1(numStrata, &tmpH);CHKERRQ(ierr);
324 ierr = PetscMalloc1(numStrata, &tmpP);CHKERRQ(ierr);
325 ierr = PetscMalloc1(numStrata, &tmpB);CHKERRQ(ierr);
326 label->numStrata = numStrata;
327 label->stratumValues = tmpV;
328 label->stratumSizes = tmpS;
329 label->ht = tmpH;
330 label->points = tmpP;
331 label->validIS = tmpB;
332 for (v = 0; v < numStrata; ++v) {
333 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr);
334 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);CHKERRQ(ierr);
335 ierr = PetscHMapISet(hmap, values[v], v);CHKERRQ(ierr);
336 tmpV[v] = values[v];
337 tmpS[v] = 0;
338 tmpH[v] = ht;
339 tmpP[v] = is;
340 tmpB[v] = PETSC_TRUE;
341 }
342 ierr = PetscObjectStateIncrease((PetscObject) label);CHKERRQ(ierr);
343 } else {
344 for (v = 0; v < numStrata; ++v) {
345 ierr = DMLabelAddStratum(label, values[v]);CHKERRQ(ierr);
346 }
347 }
348 ierr = PetscFree(values);CHKERRQ(ierr);
349 PetscFunctionReturn(0);
350 }
351
352 /*@
353 DMLabelAddStrataIS - Adds new stratum values in a DMLabel
354
355 Not collective
356
357 Input Parameter:
358 + label - The DMLabel
359 - valueIS - Index set with stratum values
360
361 Level: beginner
362
363 .seealso: DMLabelCreate(), DMLabelDestroy()
364 @*/
DMLabelAddStrataIS(DMLabel label,IS valueIS)365 PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
366 {
367 PetscInt numStrata;
368 const PetscInt *stratumValues;
369 PetscErrorCode ierr;
370
371 PetscFunctionBegin;
372 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
373 PetscValidHeaderSpecific(valueIS, IS_CLASSID, 2);
374 ierr = ISGetLocalSize(valueIS, &numStrata);CHKERRQ(ierr);
375 ierr = ISGetIndices(valueIS, &stratumValues);CHKERRQ(ierr);
376 ierr = DMLabelAddStrata(label, numStrata, stratumValues);CHKERRQ(ierr);
377 PetscFunctionReturn(0);
378 }
379
DMLabelView_Ascii(DMLabel label,PetscViewer viewer)380 static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
381 {
382 PetscInt v;
383 PetscMPIInt rank;
384 PetscErrorCode ierr;
385
386 PetscFunctionBegin;
387 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);CHKERRQ(ierr);
388 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
389 if (label) {
390 const char *name;
391
392 ierr = PetscObjectGetName((PetscObject) label, &name);CHKERRQ(ierr);
393 ierr = PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);CHKERRQ(ierr);
394 if (label->bt) {ierr = PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);CHKERRQ(ierr);}
395 for (v = 0; v < label->numStrata; ++v) {
396 const PetscInt value = label->stratumValues[v];
397 const PetscInt *points;
398 PetscInt p;
399
400 ierr = ISGetIndices(label->points[v], &points);CHKERRQ(ierr);
401 for (p = 0; p < label->stratumSizes[v]; ++p) {
402 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);CHKERRQ(ierr);
403 }
404 ierr = ISRestoreIndices(label->points[v],&points);CHKERRQ(ierr);
405 }
406 }
407 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
408 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
409 PetscFunctionReturn(0);
410 }
411
412 /*@C
413 DMLabelView - View the label
414
415 Collective on viewer
416
417 Input Parameters:
418 + label - The DMLabel
419 - viewer - The PetscViewer
420
421 Level: intermediate
422
423 .seealso: DMLabelCreate(), DMLabelDestroy()
424 @*/
DMLabelView(DMLabel label,PetscViewer viewer)425 PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
426 {
427 PetscBool iascii;
428 PetscErrorCode ierr;
429
430 PetscFunctionBegin;
431 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
432 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);CHKERRQ(ierr);}
433 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
434 if (label) {ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);}
435 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
436 if (iascii) {
437 ierr = DMLabelView_Ascii(label, viewer);CHKERRQ(ierr);
438 }
439 PetscFunctionReturn(0);
440 }
441
442 /*@
443 DMLabelReset - Destroys internal data structures in a DMLabel
444
445 Not collective
446
447 Input Parameter:
448 . label - The DMLabel
449
450 Level: beginner
451
452 .seealso: DMLabelDestroy(), DMLabelCreate()
453 @*/
DMLabelReset(DMLabel label)454 PetscErrorCode DMLabelReset(DMLabel label)
455 {
456 PetscInt v;
457 PetscErrorCode ierr;
458
459 PetscFunctionBegin;
460 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
461 for (v = 0; v < label->numStrata; ++v) {
462 ierr = PetscHSetIDestroy(&label->ht[v]);CHKERRQ(ierr);
463 ierr = ISDestroy(&label->points[v]);CHKERRQ(ierr);
464 }
465 label->numStrata = 0;
466 ierr = PetscFree(label->stratumValues);CHKERRQ(ierr);
467 ierr = PetscFree(label->stratumSizes);CHKERRQ(ierr);
468 ierr = PetscFree(label->ht);CHKERRQ(ierr);
469 ierr = PetscFree(label->points);CHKERRQ(ierr);
470 ierr = PetscFree(label->validIS);CHKERRQ(ierr);
471 ierr = PetscHMapIReset(label->hmap);CHKERRQ(ierr);
472 label->pStart = -1;
473 label->pEnd = -1;
474 ierr = PetscBTDestroy(&label->bt);CHKERRQ(ierr);
475 PetscFunctionReturn(0);
476 }
477
478 /*@
479 DMLabelDestroy - Destroys a DMLabel
480
481 Collective on label
482
483 Input Parameter:
484 . label - The DMLabel
485
486 Level: beginner
487
488 .seealso: DMLabelReset(), DMLabelCreate()
489 @*/
DMLabelDestroy(DMLabel * label)490 PetscErrorCode DMLabelDestroy(DMLabel *label)
491 {
492 PetscErrorCode ierr;
493
494 PetscFunctionBegin;
495 if (!*label) PetscFunctionReturn(0);
496 PetscValidHeaderSpecific((*label),DMLABEL_CLASSID,1);
497 if (--((PetscObject)(*label))->refct > 0) {*label = NULL; PetscFunctionReturn(0);}
498 ierr = DMLabelReset(*label);CHKERRQ(ierr);
499 ierr = PetscHMapIDestroy(&(*label)->hmap);CHKERRQ(ierr);
500 ierr = PetscHeaderDestroy(label);CHKERRQ(ierr);
501 PetscFunctionReturn(0);
502 }
503
504 /*@
505 DMLabelDuplicate - Duplicates a DMLabel
506
507 Collective on label
508
509 Input Parameter:
510 . label - The DMLabel
511
512 Output Parameter:
513 . labelnew - location to put new vector
514
515 Level: intermediate
516
517 .seealso: DMLabelCreate(), DMLabelDestroy()
518 @*/
DMLabelDuplicate(DMLabel label,DMLabel * labelnew)519 PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
520 {
521 const char *name;
522 PetscInt v;
523 PetscErrorCode ierr;
524
525 PetscFunctionBegin;
526 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
527 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
528 ierr = PetscObjectGetName((PetscObject) label, &name);CHKERRQ(ierr);
529 ierr = DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);CHKERRQ(ierr);
530
531 (*labelnew)->numStrata = label->numStrata;
532 (*labelnew)->defaultValue = label->defaultValue;
533 ierr = PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);CHKERRQ(ierr);
534 ierr = PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);CHKERRQ(ierr);
535 ierr = PetscMalloc1(label->numStrata, &(*labelnew)->ht);CHKERRQ(ierr);
536 ierr = PetscMalloc1(label->numStrata, &(*labelnew)->points);CHKERRQ(ierr);
537 ierr = PetscMalloc1(label->numStrata, &(*labelnew)->validIS);CHKERRQ(ierr);
538 for (v = 0; v < label->numStrata; ++v) {
539 ierr = PetscHSetICreate(&(*labelnew)->ht[v]);CHKERRQ(ierr);
540 (*labelnew)->stratumValues[v] = label->stratumValues[v];
541 (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
542 ierr = PetscObjectReference((PetscObject) (label->points[v]));CHKERRQ(ierr);
543 (*labelnew)->points[v] = label->points[v];
544 (*labelnew)->validIS[v] = PETSC_TRUE;
545 }
546 ierr = PetscHMapIDestroy(&(*labelnew)->hmap);CHKERRQ(ierr);
547 ierr = PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);CHKERRQ(ierr);
548 (*labelnew)->pStart = -1;
549 (*labelnew)->pEnd = -1;
550 (*labelnew)->bt = NULL;
551 PetscFunctionReturn(0);
552 }
553
554 /*@
555 DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
556
557 Not collective
558
559 Input Parameter:
560 . label - The DMLabel
561
562 Level: intermediate
563
564 .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
565 @*/
DMLabelComputeIndex(DMLabel label)566 PetscErrorCode DMLabelComputeIndex(DMLabel label)
567 {
568 PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
569 PetscErrorCode ierr;
570
571 PetscFunctionBegin;
572 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
573 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
574 for (v = 0; v < label->numStrata; ++v) {
575 const PetscInt *points;
576 PetscInt i;
577
578 ierr = ISGetIndices(label->points[v], &points);CHKERRQ(ierr);
579 for (i = 0; i < label->stratumSizes[v]; ++i) {
580 const PetscInt point = points[i];
581
582 pStart = PetscMin(point, pStart);
583 pEnd = PetscMax(point+1, pEnd);
584 }
585 ierr = ISRestoreIndices(label->points[v], &points);CHKERRQ(ierr);
586 }
587 label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
588 label->pEnd = pEnd;
589 ierr = DMLabelCreateIndex(label, label->pStart, label->pEnd);CHKERRQ(ierr);
590 PetscFunctionReturn(0);
591 }
592
593 /*@
594 DMLabelCreateIndex - Create an index structure for membership determination
595
596 Not collective
597
598 Input Parameters:
599 + label - The DMLabel
600 . pStart - The smallest point
601 - pEnd - The largest point + 1
602
603 Level: intermediate
604
605 .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
606 @*/
DMLabelCreateIndex(DMLabel label,PetscInt pStart,PetscInt pEnd)607 PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
608 {
609 PetscInt v;
610 PetscErrorCode ierr;
611
612 PetscFunctionBegin;
613 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
614 ierr = DMLabelDestroyIndex(label);CHKERRQ(ierr);
615 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
616 label->pStart = pStart;
617 label->pEnd = pEnd;
618 /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
619 ierr = PetscBTCreate(pEnd - pStart, &label->bt);CHKERRQ(ierr);
620 for (v = 0; v < label->numStrata; ++v) {
621 const PetscInt *points;
622 PetscInt i;
623
624 ierr = ISGetIndices(label->points[v], &points);CHKERRQ(ierr);
625 for (i = 0; i < label->stratumSizes[v]; ++i) {
626 const PetscInt point = points[i];
627
628 if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
629 ierr = PetscBTSet(label->bt, point - pStart);CHKERRQ(ierr);
630 }
631 ierr = ISRestoreIndices(label->points[v], &points);CHKERRQ(ierr);
632 }
633 PetscFunctionReturn(0);
634 }
635
636 /*@
637 DMLabelDestroyIndex - Destroy the index structure
638
639 Not collective
640
641 Input Parameter:
642 . label - the DMLabel
643
644 Level: intermediate
645
646 .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
647 @*/
DMLabelDestroyIndex(DMLabel label)648 PetscErrorCode DMLabelDestroyIndex(DMLabel label)
649 {
650 PetscErrorCode ierr;
651
652 PetscFunctionBegin;
653 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
654 label->pStart = -1;
655 label->pEnd = -1;
656 ierr = PetscBTDestroy(&label->bt);CHKERRQ(ierr);
657 PetscFunctionReturn(0);
658 }
659
660 /*@
661 DMLabelGetBounds - Return the smallest and largest point in the label
662
663 Not collective
664
665 Input Parameter:
666 . label - the DMLabel
667
668 Output Parameters:
669 + pStart - The smallest point
670 - pEnd - The largest point + 1
671
672 Note: This will compute an index for the label if one does not exist.
673
674 Level: intermediate
675
676 .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
677 @*/
DMLabelGetBounds(DMLabel label,PetscInt * pStart,PetscInt * pEnd)678 PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
679 {
680 PetscErrorCode ierr;
681
682 PetscFunctionBegin;
683 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
684 if ((label->pStart == -1) && (label->pEnd == -1)) {ierr = DMLabelComputeIndex(label);CHKERRQ(ierr);}
685 if (pStart) {
686 PetscValidIntPointer(pStart, 2);
687 *pStart = label->pStart;
688 }
689 if (pEnd) {
690 PetscValidIntPointer(pEnd, 3);
691 *pEnd = label->pEnd;
692 }
693 PetscFunctionReturn(0);
694 }
695
696 /*@
697 DMLabelHasValue - Determine whether a label assigns the value to any point
698
699 Not collective
700
701 Input Parameters:
702 + label - the DMLabel
703 - value - the value
704
705 Output Parameter:
706 . contains - Flag indicating whether the label maps this value to any point
707
708 Level: developer
709
710 .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
711 @*/
DMLabelHasValue(DMLabel label,PetscInt value,PetscBool * contains)712 PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
713 {
714 PetscInt v;
715 PetscErrorCode ierr;
716
717 PetscFunctionBegin;
718 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
719 PetscValidBoolPointer(contains, 3);
720 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
721 *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
722 PetscFunctionReturn(0);
723 }
724
725 /*@
726 DMLabelHasPoint - Determine whether a label assigns a value to a point
727
728 Not collective
729
730 Input Parameters:
731 + label - the DMLabel
732 - point - the point
733
734 Output Parameter:
735 . contains - Flag indicating whether the label maps this point to a value
736
737 Note: The user must call DMLabelCreateIndex() before this function.
738
739 Level: developer
740
741 .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
742 @*/
DMLabelHasPoint(DMLabel label,PetscInt point,PetscBool * contains)743 PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
744 {
745 PetscErrorCode ierr;
746
747 PetscFunctionBeginHot;
748 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
749 PetscValidBoolPointer(contains, 3);
750 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
751 if (PetscDefined(USE_DEBUG)) {
752 if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
753 if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
754 }
755 *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
756 PetscFunctionReturn(0);
757 }
758
759 /*@
760 DMLabelStratumHasPoint - Return true if the stratum contains a point
761
762 Not collective
763
764 Input Parameters:
765 + label - the DMLabel
766 . value - the stratum value
767 - point - the point
768
769 Output Parameter:
770 . contains - true if the stratum contains the point
771
772 Level: intermediate
773
774 .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
775 @*/
DMLabelStratumHasPoint(DMLabel label,PetscInt value,PetscInt point,PetscBool * contains)776 PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
777 {
778 PetscInt v;
779 PetscErrorCode ierr;
780
781 PetscFunctionBeginHot;
782 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
783 PetscValidBoolPointer(contains, 4);
784 *contains = PETSC_FALSE;
785 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
786 if (v < 0) PetscFunctionReturn(0);
787
788 if (label->validIS[v]) {
789 PetscInt i;
790
791 ierr = ISLocate(label->points[v], point, &i);CHKERRQ(ierr);
792 if (i >= 0) *contains = PETSC_TRUE;
793 } else {
794 PetscBool has;
795
796 ierr = PetscHSetIHas(label->ht[v], point, &has);CHKERRQ(ierr);
797 if (has) *contains = PETSC_TRUE;
798 }
799 PetscFunctionReturn(0);
800 }
801
802 /*@
803 DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
804 When a label is created, it is initialized to -1.
805
806 Not collective
807
808 Input parameter:
809 . label - a DMLabel object
810
811 Output parameter:
812 . defaultValue - the default value
813
814 Level: beginner
815
816 .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
817 @*/
DMLabelGetDefaultValue(DMLabel label,PetscInt * defaultValue)818 PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
819 {
820 PetscFunctionBegin;
821 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
822 *defaultValue = label->defaultValue;
823 PetscFunctionReturn(0);
824 }
825
826 /*@
827 DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
828 When a label is created, it is initialized to -1.
829
830 Not collective
831
832 Input parameter:
833 . label - a DMLabel object
834
835 Output parameter:
836 . defaultValue - the default value
837
838 Level: beginner
839
840 .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
841 @*/
DMLabelSetDefaultValue(DMLabel label,PetscInt defaultValue)842 PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
843 {
844 PetscFunctionBegin;
845 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
846 label->defaultValue = defaultValue;
847 PetscFunctionReturn(0);
848 }
849
850 /*@
851 DMLabelGetValue - Return the value a label assigns to a point, or the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue())
852
853 Not collective
854
855 Input Parameters:
856 + label - the DMLabel
857 - point - the point
858
859 Output Parameter:
860 . value - The point value, or the default value (-1 by default)
861
862 Level: intermediate
863
864 .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
865 @*/
DMLabelGetValue(DMLabel label,PetscInt point,PetscInt * value)866 PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
867 {
868 PetscInt v;
869 PetscErrorCode ierr;
870
871 PetscFunctionBeginHot;
872 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
873 PetscValidPointer(value, 3);
874 *value = label->defaultValue;
875 for (v = 0; v < label->numStrata; ++v) {
876 if (label->validIS[v]) {
877 PetscInt i;
878
879 ierr = ISLocate(label->points[v], point, &i);CHKERRQ(ierr);
880 if (i >= 0) {
881 *value = label->stratumValues[v];
882 break;
883 }
884 } else {
885 PetscBool has;
886
887 ierr = PetscHSetIHas(label->ht[v], point, &has);CHKERRQ(ierr);
888 if (has) {
889 *value = label->stratumValues[v];
890 break;
891 }
892 }
893 }
894 PetscFunctionReturn(0);
895 }
896
897 /*@
898 DMLabelSetValue - Set the value a label assigns to a point. If the value is the same as the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue() to something different), then this function will do nothing.
899
900 Not collective
901
902 Input Parameters:
903 + label - the DMLabel
904 . point - the point
905 - value - The point value
906
907 Level: intermediate
908
909 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
910 @*/
DMLabelSetValue(DMLabel label,PetscInt point,PetscInt value)911 PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
912 {
913 PetscInt v;
914 PetscErrorCode ierr;
915
916 PetscFunctionBegin;
917 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
918 /* Find label value, add new entry if needed */
919 if (value == label->defaultValue) PetscFunctionReturn(0);
920 ierr = DMLabelLookupAddStratum(label, value, &v);CHKERRQ(ierr);
921 /* Set key */
922 ierr = DMLabelMakeInvalid_Private(label, v);CHKERRQ(ierr);
923 ierr = PetscHSetIAdd(label->ht[v], point);CHKERRQ(ierr);
924 PetscFunctionReturn(0);
925 }
926
927 /*@
928 DMLabelClearValue - Clear the value a label assigns to a point
929
930 Not collective
931
932 Input Parameters:
933 + label - the DMLabel
934 . point - the point
935 - value - The point value
936
937 Level: intermediate
938
939 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
940 @*/
DMLabelClearValue(DMLabel label,PetscInt point,PetscInt value)941 PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
942 {
943 PetscInt v;
944 PetscErrorCode ierr;
945
946 PetscFunctionBegin;
947 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
948 /* Find label value */
949 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
950 if (v < 0) PetscFunctionReturn(0);
951
952 if (label->bt) {
953 if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
954 ierr = PetscBTClear(label->bt, point - label->pStart);CHKERRQ(ierr);
955 }
956
957 /* Delete key */
958 ierr = DMLabelMakeInvalid_Private(label, v);CHKERRQ(ierr);
959 ierr = PetscHSetIDel(label->ht[v], point);CHKERRQ(ierr);
960 PetscFunctionReturn(0);
961 }
962
963 /*@
964 DMLabelInsertIS - Set all points in the IS to a value
965
966 Not collective
967
968 Input Parameters:
969 + label - the DMLabel
970 . is - the point IS
971 - value - The point value
972
973 Level: intermediate
974
975 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
976 @*/
DMLabelInsertIS(DMLabel label,IS is,PetscInt value)977 PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
978 {
979 PetscInt v, n, p;
980 const PetscInt *points;
981 PetscErrorCode ierr;
982
983 PetscFunctionBegin;
984 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
985 PetscValidHeaderSpecific(is, IS_CLASSID, 2);
986 /* Find label value, add new entry if needed */
987 if (value == label->defaultValue) PetscFunctionReturn(0);
988 ierr = DMLabelLookupAddStratum(label, value, &v);CHKERRQ(ierr);
989 /* Set keys */
990 ierr = DMLabelMakeInvalid_Private(label, v);CHKERRQ(ierr);
991 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr);
992 ierr = ISGetIndices(is, &points);CHKERRQ(ierr);
993 for (p = 0; p < n; ++p) {ierr = PetscHSetIAdd(label->ht[v], points[p]);CHKERRQ(ierr);}
994 ierr = ISRestoreIndices(is, &points);CHKERRQ(ierr);
995 PetscFunctionReturn(0);
996 }
997
998 /*@
999 DMLabelGetNumValues - Get the number of values that the DMLabel takes
1000
1001 Not collective
1002
1003 Input Parameter:
1004 . label - the DMLabel
1005
1006 Output Paramater:
1007 . numValues - the number of values
1008
1009 Level: intermediate
1010
1011 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1012 @*/
DMLabelGetNumValues(DMLabel label,PetscInt * numValues)1013 PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1014 {
1015 PetscFunctionBegin;
1016 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1017 PetscValidPointer(numValues, 2);
1018 *numValues = label->numStrata;
1019 PetscFunctionReturn(0);
1020 }
1021
1022 /*@
1023 DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1024
1025 Not collective
1026
1027 Input Parameter:
1028 . label - the DMLabel
1029
1030 Output Paramater:
1031 . is - the value IS
1032
1033 Level: intermediate
1034
1035 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1036 @*/
DMLabelGetValueIS(DMLabel label,IS * values)1037 PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1038 {
1039 PetscErrorCode ierr;
1040
1041 PetscFunctionBegin;
1042 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1043 PetscValidPointer(values, 2);
1044 ierr = ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);CHKERRQ(ierr);
1045 PetscFunctionReturn(0);
1046 }
1047
1048 /*@
1049 DMLabelHasStratum - Determine whether points exist with the given value
1050
1051 Not collective
1052
1053 Input Parameters:
1054 + label - the DMLabel
1055 - value - the stratum value
1056
1057 Output Paramater:
1058 . exists - Flag saying whether points exist
1059
1060 Level: intermediate
1061
1062 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1063 @*/
DMLabelHasStratum(DMLabel label,PetscInt value,PetscBool * exists)1064 PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1065 {
1066 PetscInt v;
1067 PetscErrorCode ierr;
1068
1069 PetscFunctionBegin;
1070 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1071 PetscValidPointer(exists, 3);
1072 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
1073 *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1074 PetscFunctionReturn(0);
1075 }
1076
1077 /*@
1078 DMLabelGetStratumSize - Get the size of a stratum
1079
1080 Not collective
1081
1082 Input Parameters:
1083 + label - the DMLabel
1084 - value - the stratum value
1085
1086 Output Paramater:
1087 . size - The number of points in the stratum
1088
1089 Level: intermediate
1090
1091 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1092 @*/
DMLabelGetStratumSize(DMLabel label,PetscInt value,PetscInt * size)1093 PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1094 {
1095 PetscInt v;
1096 PetscErrorCode ierr;
1097
1098 PetscFunctionBegin;
1099 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1100 PetscValidPointer(size, 3);
1101 *size = 0;
1102 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
1103 if (v < 0) PetscFunctionReturn(0);
1104 ierr = DMLabelMakeValid_Private(label, v);CHKERRQ(ierr);
1105 *size = label->stratumSizes[v];
1106 PetscFunctionReturn(0);
1107 }
1108
1109 /*@
1110 DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1111
1112 Not collective
1113
1114 Input Parameters:
1115 + label - the DMLabel
1116 - value - the stratum value
1117
1118 Output Paramaters:
1119 + start - the smallest point in the stratum
1120 - end - the largest point in the stratum
1121
1122 Level: intermediate
1123
1124 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1125 @*/
DMLabelGetStratumBounds(DMLabel label,PetscInt value,PetscInt * start,PetscInt * end)1126 PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1127 {
1128 PetscInt v, min, max;
1129 PetscErrorCode ierr;
1130
1131 PetscFunctionBegin;
1132 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1133 if (start) {PetscValidPointer(start, 3); *start = label->defaultValue;}
1134 if (end) {PetscValidPointer(end, 4); *end = label->defaultValue;}
1135 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
1136 if (v < 0) PetscFunctionReturn(0);
1137 ierr = DMLabelMakeValid_Private(label, v);CHKERRQ(ierr);
1138 if (label->stratumSizes[v] <= 0) PetscFunctionReturn(0);
1139 ierr = ISGetMinMax(label->points[v], &min, &max);CHKERRQ(ierr);
1140 if (start) *start = min;
1141 if (end) *end = max+1;
1142 PetscFunctionReturn(0);
1143 }
1144
1145 /*@
1146 DMLabelGetStratumIS - Get an IS with the stratum points
1147
1148 Not collective
1149
1150 Input Parameters:
1151 + label - the DMLabel
1152 - value - the stratum value
1153
1154 Output Paramater:
1155 . points - The stratum points
1156
1157 Level: intermediate
1158
1159 Notes:
1160 The output IS should be destroyed when no longer needed.
1161 Returns NULL if the stratum is empty.
1162
1163 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1164 @*/
DMLabelGetStratumIS(DMLabel label,PetscInt value,IS * points)1165 PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1166 {
1167 PetscInt v;
1168 PetscErrorCode ierr;
1169
1170 PetscFunctionBegin;
1171 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1172 PetscValidPointer(points, 3);
1173 *points = NULL;
1174 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
1175 if (v < 0) PetscFunctionReturn(0);
1176 ierr = DMLabelMakeValid_Private(label, v);CHKERRQ(ierr);
1177 ierr = PetscObjectReference((PetscObject) label->points[v]);CHKERRQ(ierr);
1178 *points = label->points[v];
1179 PetscFunctionReturn(0);
1180 }
1181
1182 /*@
1183 DMLabelSetStratumIS - Set the stratum points using an IS
1184
1185 Not collective
1186
1187 Input Parameters:
1188 + label - the DMLabel
1189 . value - the stratum value
1190 - points - The stratum points
1191
1192 Level: intermediate
1193
1194 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1195 @*/
DMLabelSetStratumIS(DMLabel label,PetscInt value,IS is)1196 PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1197 {
1198 PetscInt v;
1199 PetscErrorCode ierr;
1200
1201 PetscFunctionBegin;
1202 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1203 PetscValidHeaderSpecific(is, IS_CLASSID, 3);
1204 ierr = DMLabelLookupAddStratum(label, value, &v);CHKERRQ(ierr);
1205 if (is == label->points[v]) PetscFunctionReturn(0);
1206 ierr = DMLabelClearStratum(label, value);CHKERRQ(ierr);
1207 ierr = ISGetLocalSize(is, &(label->stratumSizes[v]));CHKERRQ(ierr);
1208 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr);
1209 ierr = ISDestroy(&(label->points[v]));CHKERRQ(ierr);
1210 label->points[v] = is;
1211 label->validIS[v] = PETSC_TRUE;
1212 ierr = PetscObjectStateIncrease((PetscObject) label);CHKERRQ(ierr);
1213 if (label->bt) {
1214 const PetscInt *points;
1215 PetscInt p;
1216
1217 ierr = ISGetIndices(is,&points);CHKERRQ(ierr);
1218 for (p = 0; p < label->stratumSizes[v]; ++p) {
1219 const PetscInt point = points[p];
1220
1221 if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
1222 ierr = PetscBTSet(label->bt, point - label->pStart);CHKERRQ(ierr);
1223 }
1224 }
1225 PetscFunctionReturn(0);
1226 }
1227
1228 /*@
1229 DMLabelClearStratum - Remove a stratum
1230
1231 Not collective
1232
1233 Input Parameters:
1234 + label - the DMLabel
1235 - value - the stratum value
1236
1237 Level: intermediate
1238
1239 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1240 @*/
DMLabelClearStratum(DMLabel label,PetscInt value)1241 PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1242 {
1243 PetscInt v;
1244 PetscErrorCode ierr;
1245
1246 PetscFunctionBegin;
1247 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1248 ierr = DMLabelLookupStratum(label, value, &v);CHKERRQ(ierr);
1249 if (v < 0) PetscFunctionReturn(0);
1250 if (label->validIS[v]) {
1251 if (label->bt) {
1252 PetscInt i;
1253 const PetscInt *points;
1254
1255 ierr = ISGetIndices(label->points[v], &points);CHKERRQ(ierr);
1256 for (i = 0; i < label->stratumSizes[v]; ++i) {
1257 const PetscInt point = points[i];
1258
1259 if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
1260 ierr = PetscBTClear(label->bt, point - label->pStart);CHKERRQ(ierr);
1261 }
1262 ierr = ISRestoreIndices(label->points[v], &points);CHKERRQ(ierr);
1263 }
1264 label->stratumSizes[v] = 0;
1265 ierr = ISDestroy(&label->points[v]);CHKERRQ(ierr);
1266 ierr = ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);CHKERRQ(ierr);
1267 ierr = PetscObjectSetName((PetscObject) label->points[v], "indices");CHKERRQ(ierr);
1268 ierr = PetscObjectStateIncrease((PetscObject) label);CHKERRQ(ierr);
1269 } else {
1270 ierr = PetscHSetIClear(label->ht[v]);CHKERRQ(ierr);
1271 }
1272 PetscFunctionReturn(0);
1273 }
1274
1275 /*@
1276 DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1277
1278 Not collective
1279
1280 Input Parameters:
1281 + label - The DMLabel
1282 . value - The label value for all points
1283 . pStart - The first point
1284 - pEnd - A point beyond all marked points
1285
1286 Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1287
1288 Level: intermediate
1289
1290 .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1291 @*/
DMLabelSetStratumBounds(DMLabel label,PetscInt value,PetscInt pStart,PetscInt pEnd)1292 PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1293 {
1294 IS pIS;
1295 PetscErrorCode ierr;
1296
1297 PetscFunctionBegin;
1298 ierr = ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);CHKERRQ(ierr);
1299 ierr = DMLabelSetStratumIS(label, value, pIS);CHKERRQ(ierr);
1300 ierr = ISDestroy(&pIS);CHKERRQ(ierr);
1301 PetscFunctionReturn(0);
1302 }
1303
1304 /*@
1305 DMLabelFilter - Remove all points outside of [start, end)
1306
1307 Not collective
1308
1309 Input Parameters:
1310 + label - the DMLabel
1311 . start - the first point kept
1312 - end - one more than the last point kept
1313
1314 Level: intermediate
1315
1316 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1317 @*/
DMLabelFilter(DMLabel label,PetscInt start,PetscInt end)1318 PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1319 {
1320 PetscInt v;
1321 PetscErrorCode ierr;
1322
1323 PetscFunctionBegin;
1324 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1325 ierr = DMLabelDestroyIndex(label);CHKERRQ(ierr);
1326 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
1327 for (v = 0; v < label->numStrata; ++v) {
1328 ierr = ISGeneralFilter(label->points[v], start, end);CHKERRQ(ierr);
1329 }
1330 ierr = DMLabelCreateIndex(label, start, end);CHKERRQ(ierr);
1331 PetscFunctionReturn(0);
1332 }
1333
1334 /*@
1335 DMLabelPermute - Create a new label with permuted points
1336
1337 Not collective
1338
1339 Input Parameters:
1340 + label - the DMLabel
1341 - permutation - the point permutation
1342
1343 Output Parameter:
1344 . labelnew - the new label containing the permuted points
1345
1346 Level: intermediate
1347
1348 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1349 @*/
DMLabelPermute(DMLabel label,IS permutation,DMLabel * labelNew)1350 PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1351 {
1352 const PetscInt *perm;
1353 PetscInt numValues, numPoints, v, q;
1354 PetscErrorCode ierr;
1355
1356 PetscFunctionBegin;
1357 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1358 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2);
1359 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
1360 ierr = DMLabelDuplicate(label, labelNew);CHKERRQ(ierr);
1361 ierr = DMLabelGetNumValues(*labelNew, &numValues);CHKERRQ(ierr);
1362 ierr = ISGetLocalSize(permutation, &numPoints);CHKERRQ(ierr);
1363 ierr = ISGetIndices(permutation, &perm);CHKERRQ(ierr);
1364 for (v = 0; v < numValues; ++v) {
1365 const PetscInt size = (*labelNew)->stratumSizes[v];
1366 const PetscInt *points;
1367 PetscInt *pointsNew;
1368
1369 ierr = ISGetIndices((*labelNew)->points[v],&points);CHKERRQ(ierr);
1370 ierr = PetscMalloc1(size,&pointsNew);CHKERRQ(ierr);
1371 for (q = 0; q < size; ++q) {
1372 const PetscInt point = points[q];
1373
1374 if ((point < 0) || (point >= numPoints)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [0, %D) for the remapping", point, numPoints);
1375 pointsNew[q] = perm[point];
1376 }
1377 ierr = ISRestoreIndices((*labelNew)->points[v],&points);CHKERRQ(ierr);
1378 ierr = PetscSortInt(size, pointsNew);CHKERRQ(ierr);
1379 ierr = ISDestroy(&((*labelNew)->points[v]));CHKERRQ(ierr);
1380 if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1381 ierr = ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));CHKERRQ(ierr);
1382 ierr = PetscFree(pointsNew);CHKERRQ(ierr);
1383 } else {
1384 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));CHKERRQ(ierr);
1385 }
1386 ierr = PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");CHKERRQ(ierr);
1387 }
1388 ierr = ISRestoreIndices(permutation, &perm);CHKERRQ(ierr);
1389 if (label->bt) {
1390 ierr = PetscBTDestroy(&label->bt);CHKERRQ(ierr);
1391 ierr = DMLabelCreateIndex(label, label->pStart, label->pEnd);CHKERRQ(ierr);
1392 }
1393 PetscFunctionReturn(0);
1394 }
1395
DMLabelDistribute_Internal(DMLabel label,PetscSF sf,PetscSection * leafSection,PetscInt ** leafStrata)1396 PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1397 {
1398 MPI_Comm comm;
1399 PetscInt s, l, nroots, nleaves, dof, offset, size;
1400 PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1401 PetscSection rootSection;
1402 PetscSF labelSF;
1403 PetscErrorCode ierr;
1404
1405 PetscFunctionBegin;
1406 if (label) {ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);}
1407 ierr = PetscObjectGetComm((PetscObject)sf, &comm);CHKERRQ(ierr);
1408 /* Build a section of stratum values per point, generate the according SF
1409 and distribute point-wise stratum values to leaves. */
1410 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr);
1411 ierr = PetscSectionCreate(comm, &rootSection);CHKERRQ(ierr);
1412 ierr = PetscSectionSetChart(rootSection, 0, nroots);CHKERRQ(ierr);
1413 if (label) {
1414 for (s = 0; s < label->numStrata; ++s) {
1415 const PetscInt *points;
1416
1417 ierr = ISGetIndices(label->points[s], &points);CHKERRQ(ierr);
1418 for (l = 0; l < label->stratumSizes[s]; l++) {
1419 ierr = PetscSectionGetDof(rootSection, points[l], &dof);CHKERRQ(ierr);
1420 ierr = PetscSectionSetDof(rootSection, points[l], dof+1);CHKERRQ(ierr);
1421 }
1422 ierr = ISRestoreIndices(label->points[s], &points);CHKERRQ(ierr);
1423 }
1424 }
1425 ierr = PetscSectionSetUp(rootSection);CHKERRQ(ierr);
1426 /* Create a point-wise array of stratum values */
1427 ierr = PetscSectionGetStorageSize(rootSection, &size);CHKERRQ(ierr);
1428 ierr = PetscMalloc1(size, &rootStrata);CHKERRQ(ierr);
1429 ierr = PetscCalloc1(nroots, &rootIdx);CHKERRQ(ierr);
1430 if (label) {
1431 for (s = 0; s < label->numStrata; ++s) {
1432 const PetscInt *points;
1433
1434 ierr = ISGetIndices(label->points[s], &points);CHKERRQ(ierr);
1435 for (l = 0; l < label->stratumSizes[s]; l++) {
1436 const PetscInt p = points[l];
1437 ierr = PetscSectionGetOffset(rootSection, p, &offset);CHKERRQ(ierr);
1438 rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1439 }
1440 ierr = ISRestoreIndices(label->points[s], &points);CHKERRQ(ierr);
1441 }
1442 }
1443 /* Build SF that maps label points to remote processes */
1444 ierr = PetscSectionCreate(comm, leafSection);CHKERRQ(ierr);
1445 ierr = PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);CHKERRQ(ierr);
1446 ierr = PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);CHKERRQ(ierr);
1447 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
1448 /* Send the strata for each point over the derived SF */
1449 ierr = PetscSectionGetStorageSize(*leafSection, &size);CHKERRQ(ierr);
1450 ierr = PetscMalloc1(size, leafStrata);CHKERRQ(ierr);
1451 ierr = PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);CHKERRQ(ierr);
1452 ierr = PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);CHKERRQ(ierr);
1453 /* Clean up */
1454 ierr = PetscFree(rootStrata);CHKERRQ(ierr);
1455 ierr = PetscFree(rootIdx);CHKERRQ(ierr);
1456 ierr = PetscSectionDestroy(&rootSection);CHKERRQ(ierr);
1457 ierr = PetscSFDestroy(&labelSF);CHKERRQ(ierr);
1458 PetscFunctionReturn(0);
1459 }
1460
1461 /*@
1462 DMLabelDistribute - Create a new label pushed forward over the PetscSF
1463
1464 Collective on sf
1465
1466 Input Parameters:
1467 + label - the DMLabel
1468 - sf - the map from old to new distribution
1469
1470 Output Parameter:
1471 . labelnew - the new redistributed label
1472
1473 Level: intermediate
1474
1475 .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1476 @*/
DMLabelDistribute(DMLabel label,PetscSF sf,DMLabel * labelNew)1477 PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1478 {
1479 MPI_Comm comm;
1480 PetscSection leafSection;
1481 PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1482 PetscInt *leafStrata, *strataIdx;
1483 PetscInt **points;
1484 const char *lname = NULL;
1485 char *name;
1486 PetscInt nameSize;
1487 PetscHSetI stratumHash;
1488 size_t len = 0;
1489 PetscMPIInt rank;
1490 PetscErrorCode ierr;
1491
1492 PetscFunctionBegin;
1493 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2);
1494 if (label) {
1495 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1496 ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
1497 }
1498 ierr = PetscObjectGetComm((PetscObject)sf, &comm);CHKERRQ(ierr);
1499 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1500 /* Bcast name */
1501 if (!rank) {
1502 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr);
1503 ierr = PetscStrlen(lname, &len);CHKERRQ(ierr);
1504 }
1505 nameSize = len;
1506 ierr = MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
1507 ierr = PetscMalloc1(nameSize+1, &name);CHKERRQ(ierr);
1508 if (!rank) {ierr = PetscArraycpy(name, lname, nameSize+1);CHKERRQ(ierr);}
1509 ierr = MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);CHKERRQ(ierr);
1510 ierr = DMLabelCreate(PETSC_COMM_SELF, name, labelNew);CHKERRQ(ierr);
1511 ierr = PetscFree(name);CHKERRQ(ierr);
1512 /* Bcast defaultValue */
1513 if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1514 ierr = MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
1515 /* Distribute stratum values over the SF and get the point mapping on the receiver */
1516 ierr = DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);CHKERRQ(ierr);
1517 /* Determine received stratum values and initialise new label*/
1518 ierr = PetscHSetICreate(&stratumHash);CHKERRQ(ierr);
1519 ierr = PetscSectionGetStorageSize(leafSection, &size);CHKERRQ(ierr);
1520 for (p = 0; p < size; ++p) {ierr = PetscHSetIAdd(stratumHash, leafStrata[p]);CHKERRQ(ierr);}
1521 ierr = PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);CHKERRQ(ierr);
1522 ierr = PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);CHKERRQ(ierr);
1523 for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1524 ierr = PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);CHKERRQ(ierr);
1525 /* Turn leafStrata into indices rather than stratum values */
1526 offset = 0;
1527 ierr = PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);CHKERRQ(ierr);
1528 ierr = PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);CHKERRQ(ierr);
1529 for (s = 0; s < (*labelNew)->numStrata; ++s) {
1530 ierr = PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);CHKERRQ(ierr);
1531 }
1532 for (p = 0; p < size; ++p) {
1533 for (s = 0; s < (*labelNew)->numStrata; ++s) {
1534 if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1535 }
1536 }
1537 /* Rebuild the point strata on the receiver */
1538 ierr = PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);CHKERRQ(ierr);
1539 ierr = PetscSectionGetChart(leafSection, &pStart, &pEnd);CHKERRQ(ierr);
1540 for (p=pStart; p<pEnd; p++) {
1541 ierr = PetscSectionGetDof(leafSection, p, &dof);CHKERRQ(ierr);
1542 ierr = PetscSectionGetOffset(leafSection, p, &offset);CHKERRQ(ierr);
1543 for (s=0; s<dof; s++) {
1544 (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1545 }
1546 }
1547 ierr = PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);CHKERRQ(ierr);
1548 ierr = PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);CHKERRQ(ierr);
1549 ierr = PetscMalloc1((*labelNew)->numStrata,&points);CHKERRQ(ierr);
1550 for (s = 0; s < (*labelNew)->numStrata; ++s) {
1551 ierr = PetscHSetICreate(&(*labelNew)->ht[s]);CHKERRQ(ierr);
1552 ierr = PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));CHKERRQ(ierr);
1553 }
1554 /* Insert points into new strata */
1555 ierr = PetscCalloc1((*labelNew)->numStrata, &strataIdx);CHKERRQ(ierr);
1556 ierr = PetscSectionGetChart(leafSection, &pStart, &pEnd);CHKERRQ(ierr);
1557 for (p=pStart; p<pEnd; p++) {
1558 ierr = PetscSectionGetDof(leafSection, p, &dof);CHKERRQ(ierr);
1559 ierr = PetscSectionGetOffset(leafSection, p, &offset);CHKERRQ(ierr);
1560 for (s=0; s<dof; s++) {
1561 stratum = leafStrata[offset+s];
1562 points[stratum][strataIdx[stratum]++] = p;
1563 }
1564 }
1565 for (s = 0; s < (*labelNew)->numStrata; s++) {
1566 ierr = ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));CHKERRQ(ierr);
1567 ierr = PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");CHKERRQ(ierr);
1568 }
1569 ierr = PetscFree(points);CHKERRQ(ierr);
1570 ierr = PetscHSetIDestroy(&stratumHash);CHKERRQ(ierr);
1571 ierr = PetscFree(leafStrata);CHKERRQ(ierr);
1572 ierr = PetscFree(strataIdx);CHKERRQ(ierr);
1573 ierr = PetscSectionDestroy(&leafSection);CHKERRQ(ierr);
1574 PetscFunctionReturn(0);
1575 }
1576
1577 /*@
1578 DMLabelGather - Gather all label values from leafs into roots
1579
1580 Collective on sf
1581
1582 Input Parameters:
1583 + label - the DMLabel
1584 - sf - the Star Forest point communication map
1585
1586 Output Parameters:
1587 . labelNew - the new DMLabel with localised leaf values
1588
1589 Level: developer
1590
1591 Note: This is the inverse operation to DMLabelDistribute.
1592
1593 .seealso: DMLabelDistribute()
1594 @*/
DMLabelGather(DMLabel label,PetscSF sf,DMLabel * labelNew)1595 PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1596 {
1597 MPI_Comm comm;
1598 PetscSection rootSection;
1599 PetscSF sfLabel;
1600 PetscSFNode *rootPoints, *leafPoints;
1601 PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1602 const PetscInt *rootDegree, *ilocal;
1603 PetscInt *rootStrata;
1604 const char *lname;
1605 char *name;
1606 PetscInt nameSize;
1607 size_t len = 0;
1608 PetscMPIInt rank, size;
1609 PetscErrorCode ierr;
1610
1611 PetscFunctionBegin;
1612 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1613 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2);
1614 ierr = PetscObjectGetComm((PetscObject)sf, &comm);CHKERRQ(ierr);
1615 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1616 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
1617 /* Bcast name */
1618 if (!rank) {
1619 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr);
1620 ierr = PetscStrlen(lname, &len);CHKERRQ(ierr);
1621 }
1622 nameSize = len;
1623 ierr = MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
1624 ierr = PetscMalloc1(nameSize+1, &name);CHKERRQ(ierr);
1625 if (!rank) {ierr = PetscArraycpy(name, lname, nameSize+1);CHKERRQ(ierr);}
1626 ierr = MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);CHKERRQ(ierr);
1627 ierr = DMLabelCreate(PETSC_COMM_SELF, name, labelNew);CHKERRQ(ierr);
1628 ierr = PetscFree(name);CHKERRQ(ierr);
1629 /* Gather rank/index pairs of leaves into local roots to build
1630 an inverse, multi-rooted SF. Note that this ignores local leaf
1631 indexing due to the use of the multiSF in PetscSFGather. */
1632 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);CHKERRQ(ierr);
1633 ierr = PetscMalloc1(nroots, &leafPoints);CHKERRQ(ierr);
1634 for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1635 for (p = 0; p < nleaves; p++) {
1636 PetscInt ilp = ilocal ? ilocal[p] : p;
1637
1638 leafPoints[ilp].index = ilp;
1639 leafPoints[ilp].rank = rank;
1640 }
1641 ierr = PetscSFComputeDegreeBegin(sf, &rootDegree);CHKERRQ(ierr);
1642 ierr = PetscSFComputeDegreeEnd(sf, &rootDegree);CHKERRQ(ierr);
1643 for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1644 ierr = PetscMalloc1(nmultiroots, &rootPoints);CHKERRQ(ierr);
1645 ierr = PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);CHKERRQ(ierr);
1646 ierr = PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);CHKERRQ(ierr);
1647 ierr = PetscSFCreate(comm,& sfLabel);CHKERRQ(ierr);
1648 ierr = PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);CHKERRQ(ierr);
1649 /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1650 ierr = DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);CHKERRQ(ierr);
1651 /* Rebuild the point strata on the receiver */
1652 for (p = 0, idx = 0; p < nroots; p++) {
1653 for (d = 0; d < rootDegree[p]; d++) {
1654 ierr = PetscSectionGetDof(rootSection, idx+d, &dof);CHKERRQ(ierr);
1655 ierr = PetscSectionGetOffset(rootSection, idx+d, &offset);CHKERRQ(ierr);
1656 for (s = 0; s < dof; s++) {ierr = DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);CHKERRQ(ierr);}
1657 }
1658 idx += rootDegree[p];
1659 }
1660 ierr = PetscFree(leafPoints);CHKERRQ(ierr);
1661 ierr = PetscFree(rootStrata);CHKERRQ(ierr);
1662 ierr = PetscSectionDestroy(&rootSection);CHKERRQ(ierr);
1663 ierr = PetscSFDestroy(&sfLabel);CHKERRQ(ierr);
1664 PetscFunctionReturn(0);
1665 }
1666
1667 /*@
1668 DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1669
1670 Not collective
1671
1672 Input Parameter:
1673 . label - the DMLabel
1674
1675 Output Parameters:
1676 + section - the section giving offsets for each stratum
1677 - is - An IS containing all the label points
1678
1679 Level: developer
1680
1681 .seealso: DMLabelDistribute()
1682 @*/
DMLabelConvertToSection(DMLabel label,PetscSection * section,IS * is)1683 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1684 {
1685 IS vIS;
1686 const PetscInt *values;
1687 PetscInt *points;
1688 PetscInt nV, vS = 0, vE = 0, v, N;
1689 PetscErrorCode ierr;
1690
1691 PetscFunctionBegin;
1692 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
1693 ierr = DMLabelGetNumValues(label, &nV);CHKERRQ(ierr);
1694 ierr = DMLabelGetValueIS(label, &vIS);CHKERRQ(ierr);
1695 ierr = ISGetIndices(vIS, &values);CHKERRQ(ierr);
1696 if (nV) {vS = values[0]; vE = values[0]+1;}
1697 for (v = 1; v < nV; ++v) {
1698 vS = PetscMin(vS, values[v]);
1699 vE = PetscMax(vE, values[v]+1);
1700 }
1701 ierr = PetscSectionCreate(PETSC_COMM_SELF, section);CHKERRQ(ierr);
1702 ierr = PetscSectionSetChart(*section, vS, vE);CHKERRQ(ierr);
1703 for (v = 0; v < nV; ++v) {
1704 PetscInt n;
1705
1706 ierr = DMLabelGetStratumSize(label, values[v], &n);CHKERRQ(ierr);
1707 ierr = PetscSectionSetDof(*section, values[v], n);CHKERRQ(ierr);
1708 }
1709 ierr = PetscSectionSetUp(*section);CHKERRQ(ierr);
1710 ierr = PetscSectionGetStorageSize(*section, &N);CHKERRQ(ierr);
1711 ierr = PetscMalloc1(N, &points);CHKERRQ(ierr);
1712 for (v = 0; v < nV; ++v) {
1713 IS is;
1714 const PetscInt *spoints;
1715 PetscInt dof, off, p;
1716
1717 ierr = PetscSectionGetDof(*section, values[v], &dof);CHKERRQ(ierr);
1718 ierr = PetscSectionGetOffset(*section, values[v], &off);CHKERRQ(ierr);
1719 ierr = DMLabelGetStratumIS(label, values[v], &is);CHKERRQ(ierr);
1720 ierr = ISGetIndices(is, &spoints);CHKERRQ(ierr);
1721 for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1722 ierr = ISRestoreIndices(is, &spoints);CHKERRQ(ierr);
1723 ierr = ISDestroy(&is);CHKERRQ(ierr);
1724 }
1725 ierr = ISRestoreIndices(vIS, &values);CHKERRQ(ierr);
1726 ierr = ISDestroy(&vIS);CHKERRQ(ierr);
1727 ierr = ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
1728 PetscFunctionReturn(0);
1729 }
1730
1731 /*@
1732 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1733 the local section and an SF describing the section point overlap.
1734
1735 Collective on sf
1736
1737 Input Parameters:
1738 + s - The PetscSection for the local field layout
1739 . sf - The SF describing parallel layout of the section points
1740 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1741 . label - The label specifying the points
1742 - labelValue - The label stratum specifying the points
1743
1744 Output Parameter:
1745 . gsection - The PetscSection for the global field layout
1746
1747 Note: This gives negative sizes and offsets to points not owned by this process
1748
1749 Level: developer
1750
1751 .seealso: PetscSectionCreate()
1752 @*/
PetscSectionCreateGlobalSectionLabel(PetscSection s,PetscSF sf,PetscBool includeConstraints,DMLabel label,PetscInt labelValue,PetscSection * gsection)1753 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1754 {
1755 PetscInt *neg = NULL, *tmpOff = NULL;
1756 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1757 PetscErrorCode ierr;
1758
1759 PetscFunctionBegin;
1760 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1761 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2);
1762 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4);
1763 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);CHKERRQ(ierr);
1764 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1765 ierr = PetscSectionSetChart(*gsection, pStart, pEnd);CHKERRQ(ierr);
1766 ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
1767 if (nroots >= 0) {
1768 if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1769 ierr = PetscCalloc1(nroots, &neg);CHKERRQ(ierr);
1770 if (nroots > pEnd-pStart) {
1771 ierr = PetscCalloc1(nroots, &tmpOff);CHKERRQ(ierr);
1772 } else {
1773 tmpOff = &(*gsection)->atlasDof[-pStart];
1774 }
1775 }
1776 /* Mark ghost points with negative dof */
1777 for (p = pStart; p < pEnd; ++p) {
1778 PetscInt value;
1779
1780 ierr = DMLabelGetValue(label, p, &value);CHKERRQ(ierr);
1781 if (value != labelValue) continue;
1782 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1783 ierr = PetscSectionSetDof(*gsection, p, dof);CHKERRQ(ierr);
1784 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
1785 if (!includeConstraints && cdof > 0) {ierr = PetscSectionSetConstraintDof(*gsection, p, cdof);CHKERRQ(ierr);}
1786 if (neg) neg[p] = -(dof+1);
1787 }
1788 ierr = PetscSectionSetUpBC(*gsection);CHKERRQ(ierr);
1789 if (nroots >= 0) {
1790 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1791 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1792 if (nroots > pEnd-pStart) {
1793 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1794 }
1795 }
1796 /* Calculate new sizes, get proccess offset, and calculate point offsets */
1797 for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1798 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1799 (*gsection)->atlasOff[p] = off;
1800 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1801 }
1802 ierr = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));CHKERRQ(ierr);
1803 globalOff -= off;
1804 for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1805 (*gsection)->atlasOff[p] += globalOff;
1806 if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1807 }
1808 /* Put in negative offsets for ghost points */
1809 if (nroots >= 0) {
1810 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1811 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1812 if (nroots > pEnd-pStart) {
1813 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1814 }
1815 }
1816 if (nroots >= 0 && nroots > pEnd-pStart) {ierr = PetscFree(tmpOff);CHKERRQ(ierr);}
1817 ierr = PetscFree(neg);CHKERRQ(ierr);
1818 PetscFunctionReturn(0);
1819 }
1820
1821 typedef struct _n_PetscSectionSym_Label
1822 {
1823 DMLabel label;
1824 PetscCopyMode *modes;
1825 PetscInt *sizes;
1826 const PetscInt ***perms;
1827 const PetscScalar ***rots;
1828 PetscInt (*minMaxOrients)[2];
1829 PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1830 } PetscSectionSym_Label;
1831
PetscSectionSymLabelReset(PetscSectionSym sym)1832 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1833 {
1834 PetscInt i, j;
1835 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1836 PetscErrorCode ierr;
1837
1838 PetscFunctionBegin;
1839 for (i = 0; i <= sl->numStrata; i++) {
1840 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1841 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1842 if (sl->perms[i]) {ierr = PetscFree(sl->perms[i][j]);CHKERRQ(ierr);}
1843 if (sl->rots[i]) {ierr = PetscFree(sl->rots[i][j]);CHKERRQ(ierr);}
1844 }
1845 if (sl->perms[i]) {
1846 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1847
1848 ierr = PetscFree(perms);CHKERRQ(ierr);
1849 }
1850 if (sl->rots[i]) {
1851 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1852
1853 ierr = PetscFree(rots);CHKERRQ(ierr);
1854 }
1855 }
1856 }
1857 ierr = PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);CHKERRQ(ierr);
1858 ierr = DMLabelDestroy(&sl->label);CHKERRQ(ierr);
1859 sl->numStrata = 0;
1860 PetscFunctionReturn(0);
1861 }
1862
PetscSectionSymDestroy_Label(PetscSectionSym sym)1863 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1864 {
1865 PetscErrorCode ierr;
1866
1867 PetscFunctionBegin;
1868 ierr = PetscSectionSymLabelReset(sym);CHKERRQ(ierr);
1869 ierr = PetscFree(sym->data);CHKERRQ(ierr);
1870 PetscFunctionReturn(0);
1871 }
1872
PetscSectionSymView_Label(PetscSectionSym sym,PetscViewer viewer)1873 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1874 {
1875 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1876 PetscBool isAscii;
1877 DMLabel label = sl->label;
1878 const char *name;
1879 PetscErrorCode ierr;
1880
1881 PetscFunctionBegin;
1882 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);CHKERRQ(ierr);
1883 if (isAscii) {
1884 PetscInt i, j, k;
1885 PetscViewerFormat format;
1886
1887 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
1888 if (label) {
1889 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
1890 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1891 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1892 ierr = DMLabelView(label, viewer);CHKERRQ(ierr);
1893 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1894 } else {
1895 ierr = PetscObjectGetName((PetscObject) sl->label, &name);CHKERRQ(ierr);
1896 ierr = PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);CHKERRQ(ierr);
1897 }
1898 } else {
1899 ierr = PetscViewerASCIIPrintf(viewer, "No label given\n");CHKERRQ(ierr);
1900 }
1901 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1902 for (i = 0; i <= sl->numStrata; i++) {
1903 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1904
1905 if (!(sl->perms[i] || sl->rots[i])) {
1906 ierr = PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);CHKERRQ(ierr);
1907 } else {
1908 ierr = PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);CHKERRQ(ierr);
1909 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1910 ierr = PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);CHKERRQ(ierr);
1911 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1912 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1913 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1914 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1915 ierr = PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);CHKERRQ(ierr);
1916 } else {
1917 PetscInt tab;
1918
1919 ierr = PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);CHKERRQ(ierr);
1920 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1921 ierr = PetscViewerASCIIGetTab(viewer,&tab);CHKERRQ(ierr);
1922 if (sl->perms[i] && sl->perms[i][j]) {
1923 ierr = PetscViewerASCIIPrintf(viewer,"Permutation:");CHKERRQ(ierr);
1924 ierr = PetscViewerASCIISetTab(viewer,0);CHKERRQ(ierr);
1925 for (k = 0; k < sl->sizes[i]; k++) {ierr = PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);CHKERRQ(ierr);}
1926 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1927 ierr = PetscViewerASCIISetTab(viewer,tab);CHKERRQ(ierr);
1928 }
1929 if (sl->rots[i] && sl->rots[i][j]) {
1930 ierr = PetscViewerASCIIPrintf(viewer,"Rotations: ");CHKERRQ(ierr);
1931 ierr = PetscViewerASCIISetTab(viewer,0);CHKERRQ(ierr);
1932 #if defined(PETSC_USE_COMPLEX)
1933 for (k = 0; k < sl->sizes[i]; k++) {ierr = PetscViewerASCIIPrintf(viewer," %+f+i*%+f",PetscRealPart(sl->rots[i][j][k]),PetscImaginaryPart(sl->rots[i][j][k]));CHKERRQ(ierr);}
1934 #else
1935 for (k = 0; k < sl->sizes[i]; k++) {ierr = PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);CHKERRQ(ierr);}
1936 #endif
1937 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1938 ierr = PetscViewerASCIISetTab(viewer,tab);CHKERRQ(ierr);
1939 }
1940 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1941 }
1942 }
1943 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1944 }
1945 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1946 }
1947 }
1948 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1949 }
1950 PetscFunctionReturn(0);
1951 }
1952
1953 /*@
1954 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1955
1956 Logically collective on sym
1957
1958 Input parameters:
1959 + sym - the section symmetries
1960 - label - the DMLabel describing the types of points
1961
1962 Level: developer:
1963
1964 .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1965 @*/
PetscSectionSymLabelSetLabel(PetscSectionSym sym,DMLabel label)1966 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1967 {
1968 PetscSectionSym_Label *sl;
1969 PetscErrorCode ierr;
1970
1971 PetscFunctionBegin;
1972 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1);
1973 sl = (PetscSectionSym_Label *) sym->data;
1974 if (sl->label && sl->label != label) {ierr = PetscSectionSymLabelReset(sym);CHKERRQ(ierr);}
1975 if (label) {
1976 sl->label = label;
1977 ierr = PetscObjectReference((PetscObject) label);CHKERRQ(ierr);
1978 ierr = DMLabelGetNumValues(label,&sl->numStrata);CHKERRQ(ierr);
1979 ierr = PetscMalloc5(sl->numStrata+1,&sl->modes,sl->numStrata+1,&sl->sizes,sl->numStrata+1,&sl->perms,sl->numStrata+1,&sl->rots,sl->numStrata+1,&sl->minMaxOrients);CHKERRQ(ierr);
1980 ierr = PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));CHKERRQ(ierr);
1981 ierr = PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));CHKERRQ(ierr);
1982 ierr = PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));CHKERRQ(ierr);
1983 ierr = PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));CHKERRQ(ierr);
1984 ierr = PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));CHKERRQ(ierr);
1985 }
1986 PetscFunctionReturn(0);
1987 }
1988
1989 /*@C
1990 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1991
1992 Logically collective on sym
1993
1994 InputParameters:
1995 + sym - the section symmetries
1996 . stratum - the stratum value in the label that we are assigning symmetries for
1997 . size - the number of dofs for points in the stratum of the label
1998 . minOrient - the smallest orientation for a point in this stratum
1999 . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2000 . mode - how sym should copy the perms and rots arrays
2001 . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2002 - rots - NULL if there are no rotations, or (maxOrient - minOrient) sets of rotations, one for each orientation. A NULL set of orientations is the identity
2003
2004 Level: developer
2005
2006 .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2007 @*/
PetscSectionSymLabelSetStratum(PetscSectionSym sym,PetscInt stratum,PetscInt size,PetscInt minOrient,PetscInt maxOrient,PetscCopyMode mode,const PetscInt ** perms,const PetscScalar ** rots)2008 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2009 {
2010 PetscSectionSym_Label *sl;
2011 const char *name;
2012 PetscInt i, j, k;
2013 PetscErrorCode ierr;
2014
2015 PetscFunctionBegin;
2016 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1);
2017 sl = (PetscSectionSym_Label *) sym->data;
2018 if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
2019 for (i = 0; i <= sl->numStrata; i++) {
2020 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2021
2022 if (stratum == value) break;
2023 }
2024 ierr = PetscObjectGetName((PetscObject) sl->label, &name);CHKERRQ(ierr);
2025 if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2026 sl->sizes[i] = size;
2027 sl->modes[i] = mode;
2028 sl->minMaxOrients[i][0] = minOrient;
2029 sl->minMaxOrients[i][1] = maxOrient;
2030 if (mode == PETSC_COPY_VALUES) {
2031 if (perms) {
2032 PetscInt **ownPerms;
2033
2034 ierr = PetscCalloc1(maxOrient - minOrient,&ownPerms);CHKERRQ(ierr);
2035 for (j = 0; j < maxOrient-minOrient; j++) {
2036 if (perms[j]) {
2037 ierr = PetscMalloc1(size,&ownPerms[j]);CHKERRQ(ierr);
2038 for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2039 }
2040 }
2041 sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2042 }
2043 if (rots) {
2044 PetscScalar **ownRots;
2045
2046 ierr = PetscCalloc1(maxOrient - minOrient,&ownRots);CHKERRQ(ierr);
2047 for (j = 0; j < maxOrient-minOrient; j++) {
2048 if (rots[j]) {
2049 ierr = PetscMalloc1(size,&ownRots[j]);CHKERRQ(ierr);
2050 for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2051 }
2052 }
2053 sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2054 }
2055 } else {
2056 sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2057 sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2058 }
2059 PetscFunctionReturn(0);
2060 }
2061
PetscSectionSymGetPoints_Label(PetscSectionSym sym,PetscSection section,PetscInt numPoints,const PetscInt * points,const PetscInt ** perms,const PetscScalar ** rots)2062 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2063 {
2064 PetscInt i, j, numStrata;
2065 PetscSectionSym_Label *sl;
2066 DMLabel label;
2067 PetscErrorCode ierr;
2068
2069 PetscFunctionBegin;
2070 sl = (PetscSectionSym_Label *) sym->data;
2071 numStrata = sl->numStrata;
2072 label = sl->label;
2073 for (i = 0; i < numPoints; i++) {
2074 PetscInt point = points[2*i];
2075 PetscInt ornt = points[2*i+1];
2076
2077 for (j = 0; j < numStrata; j++) {
2078 if (label->validIS[j]) {
2079 PetscInt k;
2080
2081 ierr = ISLocate(label->points[j],point,&k);CHKERRQ(ierr);
2082 if (k >= 0) break;
2083 } else {
2084 PetscBool has;
2085
2086 ierr = PetscHSetIHas(label->ht[j], point, &has);CHKERRQ(ierr);
2087 if (has) break;
2088 }
2089 }
2090 if ((sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) && (ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1])) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D orientation %D not in range [%D, %D) for stratum %D",point,ornt,sl->minMaxOrients[j][0],sl->minMaxOrients[j][1],j < numStrata ? label->stratumValues[j] : label->defaultValue);
2091 if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2092 if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2093 }
2094 PetscFunctionReturn(0);
2095 }
2096
PetscSectionSymCreate_Label(PetscSectionSym sym)2097 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2098 {
2099 PetscSectionSym_Label *sl;
2100 PetscErrorCode ierr;
2101
2102 PetscFunctionBegin;
2103 ierr = PetscNewLog(sym,&sl);CHKERRQ(ierr);
2104 sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2105 sym->ops->view = PetscSectionSymView_Label;
2106 sym->ops->destroy = PetscSectionSymDestroy_Label;
2107 sym->data = (void *) sl;
2108 PetscFunctionReturn(0);
2109 }
2110
2111 /*@
2112 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2113
2114 Collective
2115
2116 Input Parameters:
2117 + comm - the MPI communicator for the new symmetry
2118 - label - the label defining the strata
2119
2120 Output Parameters:
2121 . sym - the section symmetries
2122
2123 Level: developer
2124
2125 .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2126 @*/
PetscSectionSymCreateLabel(MPI_Comm comm,DMLabel label,PetscSectionSym * sym)2127 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2128 {
2129 PetscErrorCode ierr;
2130
2131 PetscFunctionBegin;
2132 ierr = DMInitializePackage();CHKERRQ(ierr);
2133 ierr = PetscSectionSymCreate(comm,sym);CHKERRQ(ierr);
2134 ierr = PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);CHKERRQ(ierr);
2135 ierr = PetscSectionSymLabelSetLabel(*sym,label);CHKERRQ(ierr);
2136 PetscFunctionReturn(0);
2137 }
2138