Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates a `PetscSection` and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
39: PetscFunctionBegin;
40: PetscAssertPointer(s, 2);
41: PetscCall(ISInitializePackage());
43: PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView));
45: (*s)->pStart = -1;
46: (*s)->pEnd = -1;
47: (*s)->perm = NULL;
48: (*s)->pointMajor = PETSC_TRUE;
49: (*s)->includesConstraints = PETSC_TRUE;
50: (*s)->atlasDof = NULL;
51: (*s)->atlasOff = NULL;
52: (*s)->bc = NULL;
53: (*s)->bcIndices = NULL;
54: (*s)->setup = PETSC_FALSE;
55: (*s)->numFields = 0;
56: (*s)->fieldNames = NULL;
57: (*s)->field = NULL;
58: (*s)->useFieldOff = PETSC_FALSE;
59: (*s)->compNames = NULL;
60: (*s)->clObj = NULL;
61: (*s)->clHash = NULL;
62: (*s)->clSection = NULL;
63: (*s)->clPoints = NULL;
64: PetscCall(PetscSectionInvalidateMaxDof_Internal(*s));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*@
69: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
71: Collective
73: Input Parameter:
74: . section - the `PetscSection`
76: Output Parameter:
77: . newSection - the copy
79: Level: intermediate
81: Developer Notes:
82: What exactly does shallow mean in this context?
84: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
85: @*/
86: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
87: {
88: PetscSectionSym sym;
89: IS perm;
90: PetscInt numFields, f, c, pStart, pEnd, p;
92: PetscFunctionBegin;
95: PetscCall(PetscSectionReset(newSection));
96: PetscCall(PetscSectionGetNumFields(section, &numFields));
97: if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
98: for (f = 0; f < numFields; ++f) {
99: const char *fieldName = NULL, *compName = NULL;
100: PetscInt numComp = 0;
102: PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
103: PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
104: PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
105: PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
106: for (c = 0; c < numComp; ++c) {
107: PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
108: PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
109: }
110: PetscCall(PetscSectionGetFieldSym(section, f, &sym));
111: PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
112: }
113: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
114: PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
115: PetscCall(PetscSectionGetPermutation(section, &perm));
116: PetscCall(PetscSectionSetPermutation(newSection, perm));
117: PetscCall(PetscSectionGetSym(section, &sym));
118: PetscCall(PetscSectionSetSym(newSection, sym));
119: for (p = pStart; p < pEnd; ++p) {
120: PetscInt dof, cdof, fcdof = 0;
122: PetscCall(PetscSectionGetDof(section, p, &dof));
123: PetscCall(PetscSectionSetDof(newSection, p, dof));
124: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
125: if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
126: for (f = 0; f < numFields; ++f) {
127: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
128: PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
129: if (cdof) {
130: PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
131: if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
132: }
133: }
134: }
135: PetscCall(PetscSectionSetUp(newSection));
136: for (p = pStart; p < pEnd; ++p) {
137: PetscInt off, cdof, fcdof = 0;
138: const PetscInt *cInd;
140: /* Must set offsets in case they do not agree with the prefix sums */
141: PetscCall(PetscSectionGetOffset(section, p, &off));
142: PetscCall(PetscSectionSetOffset(newSection, p, off));
143: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
144: if (cdof) {
145: PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
146: PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
147: for (f = 0; f < numFields; ++f) {
148: PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
149: if (fcdof) {
150: PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
151: PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
152: }
153: }
154: }
155: }
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@
160: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
162: Collective
164: Input Parameter:
165: . section - the `PetscSection`
167: Output Parameter:
168: . newSection - the copy
170: Level: beginner
172: Developer Notes:
173: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
175: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
176: @*/
177: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
178: {
179: PetscFunctionBegin;
181: PetscAssertPointer(newSection, 2);
182: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
183: PetscCall(PetscSectionCopy(section, *newSection));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*@
188: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
190: Collective
192: Input Parameter:
193: . s - the `PetscSection`
195: Options Database Key:
196: . -petscsection_point_major - `PETSC_TRUE` for point-major order
198: Level: intermediate
200: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
201: @*/
202: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
203: {
204: PetscFunctionBegin;
206: PetscObjectOptionsBegin((PetscObject)s);
207: PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
208: /* process any options handlers added with PetscObjectAddOptionsHandler() */
209: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
210: PetscOptionsEnd();
211: PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: PetscSectionCompare - Compares two sections
218: Collective
220: Input Parameters:
221: + s1 - the first `PetscSection`
222: - s2 - the second `PetscSection`
224: Output Parameter:
225: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
227: Level: intermediate
229: Note:
230: Field names are disregarded.
232: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
233: @*/
234: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
235: {
236: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
237: const PetscInt *idx1, *idx2;
238: IS perm1, perm2;
239: PetscBool flg;
240: PetscMPIInt mflg;
242: PetscFunctionBegin;
245: PetscAssertPointer(congruent, 3);
246: flg = PETSC_FALSE;
248: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
249: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
250: *congruent = PETSC_FALSE;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
255: PetscCall(PetscSectionGetChart(s2, &n1, &n2));
256: if (pStart != n1 || pEnd != n2) goto not_congruent;
258: PetscCall(PetscSectionGetPermutation(s1, &perm1));
259: PetscCall(PetscSectionGetPermutation(s2, &perm2));
260: if (perm1 && perm2) {
261: PetscCall(ISEqual(perm1, perm2, congruent));
262: if (!(*congruent)) goto not_congruent;
263: } else if (perm1 != perm2) goto not_congruent;
265: for (p = pStart; p < pEnd; ++p) {
266: PetscCall(PetscSectionGetOffset(s1, p, &n1));
267: PetscCall(PetscSectionGetOffset(s2, p, &n2));
268: if (n1 != n2) goto not_congruent;
270: PetscCall(PetscSectionGetDof(s1, p, &n1));
271: PetscCall(PetscSectionGetDof(s2, p, &n2));
272: if (n1 != n2) goto not_congruent;
274: PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
275: PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
276: if (ncdof != n2) goto not_congruent;
278: PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
279: PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
280: PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
281: if (!(*congruent)) goto not_congruent;
282: }
284: PetscCall(PetscSectionGetNumFields(s1, &nfields));
285: PetscCall(PetscSectionGetNumFields(s2, &n2));
286: if (nfields != n2) goto not_congruent;
288: for (f = 0; f < nfields; ++f) {
289: PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
290: PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
291: if (n1 != n2) goto not_congruent;
293: for (p = pStart; p < pEnd; ++p) {
294: PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
295: PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
296: if (n1 != n2) goto not_congruent;
298: PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
299: PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
300: if (n1 != n2) goto not_congruent;
302: PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
303: PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
304: if (nfcdof != n2) goto not_congruent;
306: PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
307: PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
308: PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
309: if (!(*congruent)) goto not_congruent;
310: }
311: }
313: flg = PETSC_TRUE;
314: not_congruent:
315: PetscCall(MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
322: Not Collective
324: Input Parameter:
325: . s - the `PetscSection`
327: Output Parameter:
328: . numFields - the number of fields defined, or 0 if none were defined
330: Level: intermediate
332: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
333: @*/
334: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
335: {
336: PetscFunctionBegin;
338: PetscAssertPointer(numFields, 2);
339: *numFields = s->numFields;
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@
344: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
346: Not Collective
348: Input Parameters:
349: + s - the `PetscSection`
350: - numFields - the number of fields
352: Level: intermediate
354: Notes:
355: Calling this destroys all the information in the `PetscSection` including the chart.
357: You must call `PetscSectionSetChart()` after calling this.
359: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
360: @*/
361: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
362: {
363: PetscInt f;
365: PetscFunctionBegin;
367: PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
368: PetscCall(PetscSectionReset(s));
370: s->numFields = numFields;
371: PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
372: PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
373: PetscCall(PetscMalloc1(s->numFields, &s->compNames));
374: PetscCall(PetscMalloc1(s->numFields, &s->field));
375: for (f = 0; f < s->numFields; ++f) {
376: char name[64];
378: s->numFieldComponents[f] = 1;
380: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
381: PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
382: PetscCall(PetscStrallocpy(name, (char **)&s->fieldNames[f]));
383: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
384: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
385: PetscCall(PetscStrallocpy(name, (char **)&s->compNames[f][0]));
386: }
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@C
391: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
393: Not Collective
395: Input Parameters:
396: + s - the `PetscSection`
397: - field - the field number
399: Output Parameter:
400: . fieldName - the field name
402: Level: intermediate
404: Note:
405: Will error if the field number is out of range
407: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
408: @*/
409: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
410: {
411: PetscFunctionBegin;
413: PetscAssertPointer(fieldName, 3);
414: PetscSectionCheckValidField(field, s->numFields);
415: *fieldName = s->fieldNames[field];
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: /*@C
420: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
422: Not Collective
424: Input Parameters:
425: + s - the `PetscSection`
426: . field - the field number
427: - fieldName - the field name
429: Level: intermediate
431: Note:
432: Will error if the field number is out of range
434: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
435: @*/
436: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
437: {
438: PetscFunctionBegin;
440: if (fieldName) PetscAssertPointer(fieldName, 3);
441: PetscSectionCheckValidField(field, s->numFields);
442: PetscCall(PetscFree(s->fieldNames[field]));
443: PetscCall(PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*@C
448: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
450: Not Collective
452: Input Parameters:
453: + s - the `PetscSection`
454: . field - the field number
455: - comp - the component number
457: Output Parameter:
458: . compName - the component name
460: Level: intermediate
462: Note:
463: Will error if the field or component number do not exist
465: Developer Notes:
466: The function name should have Field in it since they are field components.
468: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
469: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
470: @*/
471: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
472: {
473: PetscFunctionBegin;
475: PetscAssertPointer(compName, 4);
476: PetscSectionCheckValidField(field, s->numFields);
477: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
478: *compName = s->compNames[field][comp];
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@C
483: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
485: Not Collective
487: Input Parameters:
488: + s - the `PetscSection`
489: . field - the field number
490: . comp - the component number
491: - compName - the component name
493: Level: advanced
495: Note:
496: Will error if the field or component number do not exist
498: Developer Notes:
499: The function name should have Field in it since they are field components.
501: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
502: `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
503: @*/
504: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
505: {
506: PetscFunctionBegin;
508: if (compName) PetscAssertPointer(compName, 4);
509: PetscSectionCheckValidField(field, s->numFields);
510: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
511: PetscCall(PetscFree(s->compNames[field][comp]));
512: PetscCall(PetscStrallocpy(compName, (char **)&s->compNames[field][comp]));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
519: Not Collective
521: Input Parameters:
522: + s - the `PetscSection`
523: - field - the field number
525: Output Parameter:
526: . numComp - the number of field components
528: Level: advanced
530: Developer Notes:
531: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
533: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
534: `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
535: @*/
536: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
537: {
538: PetscFunctionBegin;
540: PetscAssertPointer(numComp, 3);
541: PetscSectionCheckValidField(field, s->numFields);
542: *numComp = s->numFieldComponents[field];
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
549: Not Collective
551: Input Parameters:
552: + s - the `PetscSection`
553: . field - the field number
554: - numComp - the number of field components
556: Level: advanced
558: Note:
559: This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
560: components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
561: at a point in a discretization. For example, if in three dimensions the field is velocity, it will have 3 components, u, v, and w but
562: an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.
564: The value set with this function are not needed or used in `PetscSectionSetUp()`.
566: Developer Notes:
567: This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name
569: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
570: `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
571: @*/
572: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
573: {
574: PetscInt c;
576: PetscFunctionBegin;
578: PetscSectionCheckValidField(field, s->numFields);
579: if (s->compNames) {
580: for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
581: PetscCall(PetscFree(s->compNames[field]));
582: }
584: s->numFieldComponents[field] = numComp;
585: if (numComp) {
586: PetscCall(PetscMalloc1(numComp, (char ***)&s->compNames[field]));
587: for (c = 0; c < numComp; ++c) {
588: char name[64];
590: PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
591: PetscCall(PetscStrallocpy(name, (char **)&s->compNames[field][c]));
592: }
593: }
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*@
598: PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
600: Not Collective
602: Input Parameter:
603: . s - the `PetscSection`
605: Output Parameters:
606: + pStart - the first point
607: - pEnd - one past the last point
609: Level: intermediate
611: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
612: @*/
613: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
614: {
615: PetscFunctionBegin;
617: if (pStart) *pStart = s->pStart;
618: if (pEnd) *pEnd = s->pEnd;
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /*@
623: PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
625: Not Collective
627: Input Parameters:
628: + s - the `PetscSection`
629: . pStart - the first point
630: - pEnd - one past the last point, `pStart` $ \le $ `pEnd`
632: Level: intermediate
634: Notes:
635: The charts on different MPI processes may (and often do) overlap
637: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
639: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
641: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
642: @*/
643: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
644: {
645: PetscInt f;
647: PetscFunctionBegin;
649: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
650: /* Cannot Reset() because it destroys field information */
651: s->setup = PETSC_FALSE;
652: PetscCall(PetscSectionDestroy(&s->bc));
653: PetscCall(PetscFree(s->bcIndices));
654: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
656: s->pStart = pStart;
657: s->pEnd = pEnd;
658: PetscCall(PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff));
659: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
660: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*@
665: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
667: Not Collective
669: Input Parameter:
670: . s - the `PetscSection`
672: Output Parameter:
673: . perm - The permutation as an `IS`
675: Level: intermediate
677: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
678: @*/
679: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
680: {
681: PetscFunctionBegin;
683: if (perm) {
684: PetscAssertPointer(perm, 2);
685: *perm = s->perm;
686: }
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: /*@
691: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
693: Not Collective
695: Input Parameters:
696: + s - the `PetscSection`
697: - perm - the permutation of points
699: Level: intermediate
701: Notes:
702: The permutation must be provided before `PetscSectionSetUp()`.
704: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
706: Compart to `PetscSectionPermute()`
708: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
709: @*/
710: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
711: {
712: PetscFunctionBegin;
715: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
716: if (s->perm != perm) {
717: PetscCall(ISDestroy(&s->perm));
718: if (perm) {
719: s->perm = perm;
720: PetscCall(PetscObjectReference((PetscObject)s->perm));
721: }
722: }
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@
727: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
729: Not Collective
731: Input Parameter:
732: . s - the `PetscSection`
734: Output Parameter:
735: . pm - the flag for point major ordering
737: Level: intermediate
739: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
740: @*/
741: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
742: {
743: PetscFunctionBegin;
745: PetscAssertPointer(pm, 2);
746: *pm = s->pointMajor;
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: /*@
751: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
753: Not Collective
755: Input Parameters:
756: + s - the `PetscSection`
757: - pm - the flag for point major ordering
759: Level: intermediate
761: Note:
762: Field-major order is not recommended unless you are managing the entire problem yourself, since many higher-level functions in PETSc depend on point-major order.
764: Point major order means the degrees of freedom are stored as follows
765: .vb
766: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
767: for each point
768: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
769: .ve
771: Field major order means the degrees of freedom are stored as follows
772: .vb
773: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
774: for each field (started with unnamed default field)
775: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
776: .ve
778: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
779: @*/
780: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
781: {
782: PetscFunctionBegin;
784: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
785: s->pointMajor = pm;
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: /*@
790: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
791: The value is set with `PetscSectionSetIncludesConstraints()`
793: Not Collective
795: Input Parameter:
796: . s - the `PetscSection`
798: Output Parameter:
799: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
801: Level: intermediate
803: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
804: @*/
805: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
806: {
807: PetscFunctionBegin;
809: PetscAssertPointer(includesConstraints, 2);
810: *includesConstraints = s->includesConstraints;
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: /*@
815: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
817: Not Collective
819: Input Parameters:
820: + s - the `PetscSection`
821: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
823: Level: intermediate
825: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
826: @*/
827: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
828: {
829: PetscFunctionBegin;
831: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
832: s->includesConstraints = includesConstraints;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /*@
837: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
839: Not Collective
841: Input Parameters:
842: + s - the `PetscSection`
843: - point - the point
845: Output Parameter:
846: . numDof - the number of dof
848: Level: intermediate
850: Notes:
851: In a global section, this size will be negative for points not owned by this process.
853: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
855: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
856: @*/
857: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
858: {
859: PetscFunctionBeginHot;
861: PetscAssertPointer(numDof, 3);
862: if (PetscDefined(USE_DEBUG)) PetscCheck(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
863: *numDof = s->atlasDof[point - s->pStart];
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
870: Not Collective
872: Input Parameters:
873: + s - the `PetscSection`
874: . point - the point
875: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
877: Level: intermediate
879: Note:
880: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
882: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
883: @*/
884: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
885: {
886: PetscFunctionBegin;
888: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
889: s->atlasDof[point - s->pStart] = numDof;
890: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: /*@
895: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
897: Not Collective
899: Input Parameters:
900: + s - the `PetscSection`
901: . point - the point
902: - numDof - the number of additional dof
904: Level: intermediate
906: Note:
907: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
909: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
910: @*/
911: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
912: {
913: PetscFunctionBeginHot;
915: if (PetscDefined(USE_DEBUG)) PetscCheck(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
916: s->atlasDof[point - s->pStart] += numDof;
917: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: /*@
922: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
924: Not Collective
926: Input Parameters:
927: + s - the `PetscSection`
928: . point - the point
929: - field - the field
931: Output Parameter:
932: . numDof - the number of dof
934: Level: intermediate
936: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
937: @*/
938: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
939: {
940: PetscFunctionBegin;
942: PetscAssertPointer(numDof, 4);
943: PetscSectionCheckValidField(field, s->numFields);
944: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: /*@
949: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
951: Not Collective
953: Input Parameters:
954: + s - the `PetscSection`
955: . point - the point
956: . field - the field
957: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
959: Level: intermediate
961: Note:
962: When setting the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
963: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
965: This is equivalent to
966: .vb
967: PetscSection fs;
968: PetscSectionGetField(s,field,&fs)
969: PetscSectionSetDof(fs,numDof)
970: .ve
972: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
973: @*/
974: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
975: {
976: PetscFunctionBegin;
978: PetscSectionCheckValidField(field, s->numFields);
979: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
986: Not Collective
988: Input Parameters:
989: + s - the `PetscSection`
990: . point - the point
991: . field - the field
992: - numDof - the number of dof
994: Level: intermediate
996: Notes:
997: When adding to the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
998: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1000: This is equivalent to
1001: .vb
1002: PetscSection fs;
1003: PetscSectionGetField(s,field,&fs)
1004: PetscSectionAddDof(fs,numDof)
1005: .ve
1007: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1008: @*/
1009: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1010: {
1011: PetscFunctionBegin;
1013: PetscSectionCheckValidField(field, s->numFields);
1014: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@
1019: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1021: Not Collective
1023: Input Parameters:
1024: + s - the `PetscSection`
1025: - point - the point
1027: Output Parameter:
1028: . numDof - the number of dof which are fixed by constraints
1030: Level: intermediate
1032: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1033: @*/
1034: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1035: {
1036: PetscFunctionBegin;
1038: PetscAssertPointer(numDof, 3);
1039: if (s->bc) {
1040: PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1041: } else *numDof = 0;
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@
1046: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
1048: Not Collective
1050: Input Parameters:
1051: + s - the `PetscSection`
1052: . point - the point
1053: - numDof - the number of dof which are fixed by constraints
1055: Level: intermediate
1057: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1058: @*/
1059: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1060: {
1061: PetscFunctionBegin;
1063: if (numDof) {
1064: PetscCall(PetscSectionCheckConstraints_Private(s));
1065: PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1066: }
1067: PetscFunctionReturn(PETSC_SUCCESS);
1068: }
1070: /*@
1071: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1073: Not Collective
1075: Input Parameters:
1076: + s - the `PetscSection`
1077: . point - the point
1078: - numDof - the number of additional dof which are fixed by constraints
1080: Level: intermediate
1082: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1083: @*/
1084: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1085: {
1086: PetscFunctionBegin;
1088: if (numDof) {
1089: PetscCall(PetscSectionCheckConstraints_Private(s));
1090: PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1091: }
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1098: Not Collective
1100: Input Parameters:
1101: + s - the `PetscSection`
1102: . point - the point
1103: - field - the field
1105: Output Parameter:
1106: . numDof - the number of dof which are fixed by constraints
1108: Level: intermediate
1110: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1111: @*/
1112: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1113: {
1114: PetscFunctionBegin;
1116: PetscAssertPointer(numDof, 4);
1117: PetscSectionCheckValidField(field, s->numFields);
1118: PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: /*@
1123: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1125: Not Collective
1127: Input Parameters:
1128: + s - the `PetscSection`
1129: . point - the point
1130: . field - the field
1131: - numDof - the number of dof which are fixed by constraints
1133: Level: intermediate
1135: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1136: @*/
1137: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1138: {
1139: PetscFunctionBegin;
1141: PetscSectionCheckValidField(field, s->numFields);
1142: PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1143: PetscFunctionReturn(PETSC_SUCCESS);
1144: }
1146: /*@
1147: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1149: Not Collective
1151: Input Parameters:
1152: + s - the `PetscSection`
1153: . point - the point
1154: . field - the field
1155: - numDof - the number of additional dof which are fixed by constraints
1157: Level: intermediate
1159: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1160: @*/
1161: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1162: {
1163: PetscFunctionBegin;
1165: PetscSectionCheckValidField(field, s->numFields);
1166: PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: /*@
1171: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1173: Not Collective
1175: Input Parameter:
1176: . s - the `PetscSection`
1178: Level: advanced
1180: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1181: @*/
1182: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1183: {
1184: PetscFunctionBegin;
1186: if (s->bc) {
1187: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1189: PetscCall(PetscSectionSetUp(s->bc));
1190: PetscCall(PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices));
1191: }
1192: PetscFunctionReturn(PETSC_SUCCESS);
1193: }
1195: /*@
1196: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`
1198: Not Collective
1200: Input Parameter:
1201: . s - the `PetscSection`
1203: Level: intermediate
1205: Notes:
1206: If used, `PetscSectionSetPermutation()` must be called before this routine.
1208: `PetscSectionSetPointMajor()`, cannot be called after this routine.
1210: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1211: @*/
1212: PetscErrorCode PetscSectionSetUp(PetscSection s)
1213: {
1214: const PetscInt *pind = NULL;
1215: PetscInt offset = 0, foff, p, f;
1217: PetscFunctionBegin;
1219: if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1220: s->setup = PETSC_TRUE;
1221: /* Set offsets and field offsets for all points */
1222: /* Assume that all fields have the same chart */
1223: PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1224: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1225: if (s->pointMajor) {
1226: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1227: const PetscInt q = pind ? pind[p] : p;
1229: /* Set point offset */
1230: s->atlasOff[q] = offset;
1231: offset += s->atlasDof[q];
1232: /* Set field offset */
1233: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1234: PetscSection sf = s->field[f];
1236: sf->atlasOff[q] = foff;
1237: foff += sf->atlasDof[q];
1238: }
1239: }
1240: } else {
1241: /* Set field offsets for all points */
1242: for (f = 0; f < s->numFields; ++f) {
1243: PetscSection sf = s->field[f];
1245: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1246: const PetscInt q = pind ? pind[p] : p;
1248: sf->atlasOff[q] = offset;
1249: offset += sf->atlasDof[q];
1250: }
1251: }
1252: /* Disable point offsets since these are unused */
1253: for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1254: }
1255: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1256: /* Setup BC sections */
1257: PetscCall(PetscSectionSetUpBC(s));
1258: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1259: PetscFunctionReturn(PETSC_SUCCESS);
1260: }
1262: /*@
1263: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`
1265: Not Collective
1267: Input Parameter:
1268: . s - the `PetscSection`
1270: Output Parameter:
1271: . maxDof - the maximum dof
1273: Level: intermediate
1275: Notes:
1276: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1278: This is the maximum over all points of the sum of the number of dof in the unnamed default field plus all named fields. This is equivalent to
1279: the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process
1281: Developer Notes:
1282: The returned number is calculated lazily and stashed.
1284: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1286: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1288: It should also be called every time `atlasDof` is modified directly.
1290: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1291: @*/
1292: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1293: {
1294: PetscInt p;
1296: PetscFunctionBegin;
1298: PetscAssertPointer(maxDof, 2);
1299: if (s->maxDof == PETSC_MIN_INT) {
1300: s->maxDof = 0;
1301: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1302: }
1303: *maxDof = s->maxDof;
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: /*@
1308: PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`
1310: Not Collective
1312: Input Parameter:
1313: . s - the `PetscSection`
1315: Output Parameter:
1316: . size - the size of an array which can hold all the dofs
1318: Level: intermediate
1320: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1321: @*/
1322: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1323: {
1324: PetscInt p, n = 0;
1326: PetscFunctionBegin;
1328: PetscAssertPointer(size, 2);
1329: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1330: *size = n;
1331: PetscFunctionReturn(PETSC_SUCCESS);
1332: }
1334: /*@
1335: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`
1337: Not Collective
1339: Input Parameter:
1340: . s - the `PetscSection`
1342: Output Parameter:
1343: . size - the size of an array which can hold all unconstrained dofs
1345: Level: intermediate
1347: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1348: @*/
1349: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1350: {
1351: PetscInt p, n = 0;
1353: PetscFunctionBegin;
1355: PetscAssertPointer(size, 2);
1356: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1357: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1358: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1359: }
1360: *size = n;
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@
1365: PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1366: a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.
1368: Input Parameters:
1369: + s - The `PetscSection` for the local field layout
1370: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1371: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1372: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1374: Output Parameter:
1375: . gsection - The `PetscSection` for the global field layout
1377: Level: intermediate
1379: Notes:
1380: On each MPI process `gsection` inherits the chart of the `s` on that process.
1382: This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1383: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1385: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1386: @*/
1387: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1388: {
1389: PetscSection gs;
1390: const PetscInt *pind = NULL;
1391: PetscInt *recv = NULL, *neg = NULL;
1392: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1393: PetscInt numFields, f, numComponents;
1395: PetscFunctionBegin;
1400: PetscAssertPointer(gsection, 5);
1401: PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1402: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1403: PetscCall(PetscSectionGetNumFields(s, &numFields));
1404: if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1405: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1406: PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1407: gs->includesConstraints = includeConstraints;
1408: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1409: nlocal = nroots; /* The local/leaf space matches global/root space */
1410: /* Must allocate for all points visible to SF, which may be more than this section */
1411: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1412: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1413: PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1414: PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1415: PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1416: PetscCall(PetscArrayzero(neg, nroots));
1417: }
1418: /* Mark all local points with negative dof */
1419: for (p = pStart; p < pEnd; ++p) {
1420: PetscCall(PetscSectionGetDof(s, p, &dof));
1421: PetscCall(PetscSectionSetDof(gs, p, dof));
1422: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1423: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1424: if (neg) neg[p] = -(dof + 1);
1425: }
1426: PetscCall(PetscSectionSetUpBC(gs));
1427: if (gs->bcIndices) PetscCall(PetscArraycpy(gs->bcIndices, s->bcIndices, gs->bc->atlasOff[gs->bc->pEnd - gs->bc->pStart - 1] + gs->bc->atlasDof[gs->bc->pEnd - gs->bc->pStart - 1]));
1428: if (nroots >= 0) {
1429: PetscCall(PetscArrayzero(recv, nlocal));
1430: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1431: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1432: for (p = pStart; p < pEnd; ++p) {
1433: if (recv[p] < 0) {
1434: gs->atlasDof[p - pStart] = recv[p];
1435: PetscCall(PetscSectionGetDof(s, p, &dof));
1436: PetscCheck(-(recv[p] + 1) == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %" PetscInt_FMT " for point %" PetscInt_FMT " is not the unconstrained %" PetscInt_FMT, -(recv[p] + 1), p, dof);
1437: }
1438: }
1439: }
1440: /* Calculate new sizes, get process offset, and calculate point offsets */
1441: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1442: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1443: const PetscInt q = pind ? pind[p] : p;
1445: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1446: gs->atlasOff[q] = off;
1447: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1448: }
1449: if (!localOffsets) {
1450: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1451: globalOff -= off;
1452: }
1453: for (p = pStart, off = 0; p < pEnd; ++p) {
1454: gs->atlasOff[p - pStart] += globalOff;
1455: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1456: }
1457: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1458: /* Put in negative offsets for ghost points */
1459: if (nroots >= 0) {
1460: PetscCall(PetscArrayzero(recv, nlocal));
1461: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1462: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1463: for (p = pStart; p < pEnd; ++p) {
1464: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1465: }
1466: }
1467: PetscCall(PetscFree2(neg, recv));
1468: /* Set field dofs/offsets/constraints */
1469: for (f = 0; f < numFields; ++f) {
1470: gs->field[f]->includesConstraints = includeConstraints;
1471: PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1472: PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1473: }
1474: for (p = pStart; p < pEnd; ++p) {
1475: PetscCall(PetscSectionGetOffset(gs, p, &off));
1476: for (f = 0, foff = off; f < numFields; ++f) {
1477: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1478: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1479: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1480: PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1481: PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1482: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1483: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1484: }
1485: }
1486: for (f = 0; f < numFields; ++f) {
1487: PetscSection gfs = gs->field[f];
1489: PetscCall(PetscSectionSetUpBC(gfs));
1490: if (gfs->bcIndices) PetscCall(PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd - gfs->bc->pStart - 1] + gfs->bc->atlasDof[gfs->bc->pEnd - gfs->bc->pStart - 1]));
1491: }
1492: gs->setup = PETSC_TRUE;
1493: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1494: *gsection = gs;
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: /*@
1499: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1500: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1502: Input Parameters:
1503: + s - The `PetscSection` for the local field layout
1504: . sf - The `PetscSF` describing parallel layout of the section points
1505: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1506: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1507: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are `numExcludes` pairs and must have the same values on all MPI processes
1509: Output Parameter:
1510: . gsection - The `PetscSection` for the global field layout
1512: Level: advanced
1514: Notes:
1515: On each MPI process `gsection` inherits the chart of the `s` on that process.
1517: This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1518: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1520: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1522: Developer Notes:
1523: This is a terrible function name
1525: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1526: @*/
1527: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1528: {
1529: const PetscInt *pind = NULL;
1530: PetscInt *neg = NULL, *tmpOff = NULL;
1531: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1533: PetscFunctionBegin;
1536: PetscAssertPointer(gsection, 6);
1537: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1538: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1539: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1540: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1541: if (nroots >= 0) {
1542: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1543: PetscCall(PetscCalloc1(nroots, &neg));
1544: if (nroots > pEnd - pStart) {
1545: PetscCall(PetscCalloc1(nroots, &tmpOff));
1546: } else {
1547: tmpOff = &(*gsection)->atlasDof[-pStart];
1548: }
1549: }
1550: /* Mark ghost points with negative dof */
1551: for (p = pStart; p < pEnd; ++p) {
1552: for (e = 0; e < numExcludes; ++e) {
1553: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1554: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1555: break;
1556: }
1557: }
1558: if (e < numExcludes) continue;
1559: PetscCall(PetscSectionGetDof(s, p, &dof));
1560: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1561: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1562: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1563: if (neg) neg[p] = -(dof + 1);
1564: }
1565: PetscCall(PetscSectionSetUpBC(*gsection));
1566: if (nroots >= 0) {
1567: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1568: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1569: if (nroots > pEnd - pStart) {
1570: for (p = pStart; p < pEnd; ++p) {
1571: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1572: }
1573: }
1574: }
1575: /* Calculate new sizes, get process offset, and calculate point offsets */
1576: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1577: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1578: const PetscInt q = pind ? pind[p] : p;
1580: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1581: (*gsection)->atlasOff[q] = off;
1582: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1583: }
1584: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1585: globalOff -= off;
1586: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1587: (*gsection)->atlasOff[p] += globalOff;
1588: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1589: }
1590: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1591: /* Put in negative offsets for ghost points */
1592: if (nroots >= 0) {
1593: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1594: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1595: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1596: if (nroots > pEnd - pStart) {
1597: for (p = pStart; p < pEnd; ++p) {
1598: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1599: }
1600: }
1601: }
1602: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1603: PetscCall(PetscFree(neg));
1604: PetscFunctionReturn(PETSC_SUCCESS);
1605: }
1607: /*@
1608: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1610: Collective
1612: Input Parameters:
1613: + comm - The `MPI_Comm`
1614: - s - The `PetscSection`
1616: Output Parameter:
1617: . layout - The point layout for the data that defines the section
1619: Level: advanced
1621: Notes:
1622: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1623: degrees of freedom).
1625: This count includes constrained degrees of freedom
1627: This is usually called on the default global section.
1629: Example:
1630: .vb
1631: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1632: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1633: .ve
1635: Developer Notes:
1636: I find the names of these two functions extremely non-informative
1638: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1639: @*/
1640: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1641: {
1642: PetscInt pStart, pEnd, p, localSize = 0;
1644: PetscFunctionBegin;
1645: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1646: for (p = pStart; p < pEnd; ++p) {
1647: PetscInt dof;
1649: PetscCall(PetscSectionGetDof(s, p, &dof));
1650: if (dof >= 0) ++localSize;
1651: }
1652: PetscCall(PetscLayoutCreate(comm, layout));
1653: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1654: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1655: PetscCall(PetscLayoutSetUp(*layout));
1656: PetscFunctionReturn(PETSC_SUCCESS);
1657: }
1659: /*@
1660: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1662: Collective
1664: Input Parameters:
1665: + comm - The `MPI_Comm`
1666: - s - The `PetscSection`
1668: Output Parameter:
1669: . layout - The dof layout for the section
1671: Level: advanced
1673: Notes:
1674: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1675: including the constrained degrees of freedom
1677: This is usually called for the default global section.
1679: Example:
1680: .vb
1681: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1682: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1683: .ve
1685: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1686: @*/
1687: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1688: {
1689: PetscInt pStart, pEnd, p, localSize = 0;
1691: PetscFunctionBegin;
1693: PetscAssertPointer(layout, 3);
1694: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1695: for (p = pStart; p < pEnd; ++p) {
1696: PetscInt dof, cdof;
1698: PetscCall(PetscSectionGetDof(s, p, &dof));
1699: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1700: if (dof - cdof > 0) localSize += dof - cdof;
1701: }
1702: PetscCall(PetscLayoutCreate(comm, layout));
1703: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1704: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1705: PetscCall(PetscLayoutSetUp(*layout));
1706: PetscFunctionReturn(PETSC_SUCCESS);
1707: }
1709: /*@
1710: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1712: Not Collective
1714: Input Parameters:
1715: + s - the `PetscSection`
1716: - point - the point
1718: Output Parameter:
1719: . offset - the offset
1721: Level: intermediate
1723: Notes:
1724: In a global section, `offset` will be negative for points not owned by this process.
1726: This is for the unnamed default field in the `PetscSection` not the named fields
1728: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1730: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1731: @*/
1732: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1733: {
1734: PetscFunctionBegin;
1736: PetscAssertPointer(offset, 3);
1737: if (PetscDefined(USE_DEBUG)) PetscCheck(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
1738: *offset = s->atlasOff[point - s->pStart];
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1745: Not Collective
1747: Input Parameters:
1748: + s - the `PetscSection`
1749: . point - the point
1750: - offset - the offset, these values may be negative indicating the values are off process
1752: Level: developer
1754: Note:
1755: The user usually does not call this function, but uses `PetscSectionSetUp()`
1757: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1758: @*/
1759: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1760: {
1761: PetscFunctionBegin;
1763: PetscCheck(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
1764: s->atlasOff[point - s->pStart] = offset;
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: /*@
1769: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1771: Not Collective
1773: Input Parameters:
1774: + s - the `PetscSection`
1775: . point - the point
1776: - field - the field
1778: Output Parameter:
1779: . offset - the offset
1781: Level: intermediate
1783: Notes:
1784: In a global section, `offset` will be negative for points not owned by this process.
1786: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1788: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1789: @*/
1790: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1791: {
1792: PetscFunctionBegin;
1794: PetscAssertPointer(offset, 4);
1795: PetscSectionCheckValidField(field, s->numFields);
1796: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1797: PetscFunctionReturn(PETSC_SUCCESS);
1798: }
1800: /*@
1801: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1803: Not Collective
1805: Input Parameters:
1806: + s - the `PetscSection`
1807: . point - the point
1808: . field - the field
1809: - offset - the offset, these values may be negative indicating the values are off process
1811: Level: developer
1813: Note:
1814: The user usually does not call this function, but uses `PetscSectionSetUp()`
1816: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1817: @*/
1818: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1819: {
1820: PetscFunctionBegin;
1822: PetscSectionCheckValidField(field, s->numFields);
1823: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }
1827: /*@
1828: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1829: unnamed default field's first dof
1831: Not Collective
1833: Input Parameters:
1834: + s - the `PetscSection`
1835: . point - the point
1836: - field - the field
1838: Output Parameter:
1839: . offset - the offset
1841: Level: advanced
1843: Note:
1844: This ignores constraints
1846: Example:
1847: .vb
1848: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1849: The unnamed default field has 3 dof at `point`
1850: Field 0 has 2 dof at `point`
1851: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1852: .ve
1854: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1855: @*/
1856: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1857: {
1858: PetscInt off, foff;
1860: PetscFunctionBegin;
1862: PetscAssertPointer(offset, 4);
1863: PetscSectionCheckValidField(field, s->numFields);
1864: PetscCall(PetscSectionGetOffset(s, point, &off));
1865: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1866: *offset = foff - off;
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: /*@
1871: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1873: Not Collective
1875: Input Parameter:
1876: . s - the `PetscSection`
1878: Output Parameters:
1879: + start - the minimum offset
1880: - end - one more than the maximum offset
1882: Level: intermediate
1884: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1885: @*/
1886: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1887: {
1888: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1890: PetscFunctionBegin;
1892: if (s->atlasOff) {
1893: os = s->atlasOff[0];
1894: oe = s->atlasOff[0];
1895: }
1896: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1897: for (p = 0; p < pEnd - pStart; ++p) {
1898: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1900: if (off >= 0) {
1901: os = PetscMin(os, off);
1902: oe = PetscMax(oe, off + dof);
1903: }
1904: }
1905: if (start) *start = os;
1906: if (end) *end = oe;
1907: PetscFunctionReturn(PETSC_SUCCESS);
1908: }
1910: /*@
1911: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
1913: Collective
1915: Input Parameters:
1916: + s - the `PetscSection`
1917: . len - the number of subfields
1918: - fields - the subfield numbers
1920: Output Parameter:
1921: . subs - the subsection
1923: Level: advanced
1925: Notes:
1926: The chart of `subs` is the same as the chart of `s`
1928: This will error if a fieldnumber is out of range
1930: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1931: @*/
1932: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1933: {
1934: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1936: PetscFunctionBegin;
1937: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
1939: PetscAssertPointer(fields, 3);
1940: PetscAssertPointer(subs, 4);
1941: PetscCall(PetscSectionGetNumFields(s, &nF));
1942: PetscCheck(len <= nF, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of fields %" PetscInt_FMT, len, nF);
1943: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
1944: PetscCall(PetscSectionSetNumFields(*subs, len));
1945: for (f = 0; f < len; ++f) {
1946: const char *name = NULL;
1947: PetscInt numComp = 0;
1948: PetscSectionSym sym;
1950: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
1951: PetscCall(PetscSectionSetFieldName(*subs, f, name));
1952: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
1953: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
1954: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1955: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
1956: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
1957: }
1958: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
1959: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
1960: }
1961: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1962: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
1963: for (p = pStart; p < pEnd; ++p) {
1964: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1966: for (f = 0; f < len; ++f) {
1967: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1968: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
1969: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1970: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
1971: dof += fdof;
1972: cdof += cfdof;
1973: }
1974: PetscCall(PetscSectionSetDof(*subs, p, dof));
1975: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
1976: maxCdof = PetscMax(cdof, maxCdof);
1977: }
1978: PetscCall(PetscSectionSetUp(*subs));
1979: if (maxCdof) {
1980: PetscInt *indices;
1982: PetscCall(PetscMalloc1(maxCdof, &indices));
1983: for (p = pStart; p < pEnd; ++p) {
1984: PetscInt cdof;
1986: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
1987: if (cdof) {
1988: const PetscInt *oldIndices = NULL;
1989: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1991: for (f = 0; f < len; ++f) {
1992: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1993: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1994: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
1995: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
1996: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1997: numConst += cfdof;
1998: fOff += fdof;
1999: }
2000: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2001: }
2002: }
2003: PetscCall(PetscFree(indices));
2004: }
2005: PetscFunctionReturn(PETSC_SUCCESS);
2006: }
2008: /*@
2009: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2011: Collective
2013: Input Parameters:
2014: + s - the input sections
2015: - len - the number of input sections
2017: Output Parameter:
2018: . supers - the supersection
2020: Level: advanced
2022: Notes:
2023: The section offsets now refer to a new, larger vector.
2025: Developer Notes:
2026: Needs to explain how the sections are composed
2028: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2029: @*/
2030: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2031: {
2032: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
2034: PetscFunctionBegin;
2035: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2036: for (i = 0; i < len; ++i) {
2037: PetscInt nf, pStarti, pEndi;
2039: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2040: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2041: pStart = PetscMin(pStart, pStarti);
2042: pEnd = PetscMax(pEnd, pEndi);
2043: Nf += nf;
2044: }
2045: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2046: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2047: for (i = 0, f = 0; i < len; ++i) {
2048: PetscInt nf, fi, ci;
2050: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2051: for (fi = 0; fi < nf; ++fi, ++f) {
2052: const char *name = NULL;
2053: PetscInt numComp = 0;
2055: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2056: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2057: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2058: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2059: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2060: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2061: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2062: }
2063: }
2064: }
2065: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2066: for (p = pStart; p < pEnd; ++p) {
2067: PetscInt dof = 0, cdof = 0;
2069: for (i = 0, f = 0; i < len; ++i) {
2070: PetscInt nf, fi, pStarti, pEndi;
2071: PetscInt fdof = 0, cfdof = 0;
2073: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2074: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2075: if ((p < pStarti) || (p >= pEndi)) continue;
2076: for (fi = 0; fi < nf; ++fi, ++f) {
2077: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2078: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2079: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2080: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2081: dof += fdof;
2082: cdof += cfdof;
2083: }
2084: }
2085: PetscCall(PetscSectionSetDof(*supers, p, dof));
2086: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2087: maxCdof = PetscMax(cdof, maxCdof);
2088: }
2089: PetscCall(PetscSectionSetUp(*supers));
2090: if (maxCdof) {
2091: PetscInt *indices;
2093: PetscCall(PetscMalloc1(maxCdof, &indices));
2094: for (p = pStart; p < pEnd; ++p) {
2095: PetscInt cdof;
2097: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2098: if (cdof) {
2099: PetscInt dof, numConst = 0, fOff = 0;
2101: for (i = 0, f = 0; i < len; ++i) {
2102: const PetscInt *oldIndices = NULL;
2103: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2105: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2106: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2107: if ((p < pStarti) || (p >= pEndi)) continue;
2108: for (fi = 0; fi < nf; ++fi, ++f) {
2109: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2110: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2111: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2112: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2113: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2114: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2115: numConst += cfdof;
2116: }
2117: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2118: fOff += dof;
2119: }
2120: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2121: }
2122: }
2123: PetscCall(PetscFree(indices));
2124: }
2125: PetscFunctionReturn(PETSC_SUCCESS);
2126: }
2128: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
2129: {
2130: const PetscInt *points = NULL, *indices = NULL;
2131: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2133: PetscFunctionBegin;
2136: PetscAssertPointer(subs, 4);
2137: PetscCall(PetscSectionGetNumFields(s, &numFields));
2138: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2139: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2140: for (f = 0; f < numFields; ++f) {
2141: const char *name = NULL;
2142: PetscInt numComp = 0;
2144: PetscCall(PetscSectionGetFieldName(s, f, &name));
2145: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2146: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2147: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2148: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2149: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2150: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2151: }
2152: }
2153: /* For right now, we do not try to squeeze the subchart */
2154: if (subpointMap) {
2155: PetscCall(ISGetSize(subpointMap, &numSubpoints));
2156: PetscCall(ISGetIndices(subpointMap, &points));
2157: }
2158: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2159: if (renumberPoints) {
2160: spStart = 0;
2161: spEnd = numSubpoints;
2162: } else {
2163: PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd));
2164: ++spEnd;
2165: }
2166: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2167: for (p = pStart; p < pEnd; ++p) {
2168: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2170: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2171: if (subp < 0) continue;
2172: if (!renumberPoints) subp = p;
2173: for (f = 0; f < numFields; ++f) {
2174: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2175: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2176: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2177: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2178: }
2179: PetscCall(PetscSectionGetDof(s, p, &dof));
2180: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2181: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2182: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2183: }
2184: PetscCall(PetscSectionSetUp(*subs));
2185: /* Change offsets to original offsets */
2186: for (p = pStart; p < pEnd; ++p) {
2187: PetscInt off, foff = 0;
2189: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2190: if (subp < 0) continue;
2191: if (!renumberPoints) subp = p;
2192: for (f = 0; f < numFields; ++f) {
2193: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2194: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2195: }
2196: PetscCall(PetscSectionGetOffset(s, p, &off));
2197: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2198: }
2199: /* Copy constraint indices */
2200: for (subp = spStart; subp < spEnd; ++subp) {
2201: PetscInt cdof;
2203: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2204: if (cdof) {
2205: for (f = 0; f < numFields; ++f) {
2206: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2207: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2208: }
2209: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2210: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2211: }
2212: }
2213: if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points));
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: /*@
2218: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2220: Collective
2222: Input Parameters:
2223: + s - the `PetscSection`
2224: - subpointMap - a sorted list of points in the original mesh which are in the submesh
2226: Output Parameter:
2227: . subs - the subsection
2229: Level: advanced
2231: Notes:
2232: The points are renumbered from 0, and the section offsets now refer to a new, smaller vector. That is the chart of `subs` is `[0,sizeof(subpointmap))`
2234: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2236: Developer Notes:
2237: The use of the term Submesh is confusing and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`
2239: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2240: @*/
2241: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2242: {
2243: PetscFunctionBegin;
2244: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_TRUE, subs));
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: /*@
2249: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2251: Collective
2253: Input Parameters:
2254: + s - the `PetscSection`
2255: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2257: Output Parameter:
2258: . subs - the subsection
2260: Level: advanced
2262: Notes:
2263: The point numbers remain the same as in the larger `PetscSection`, but the section offsets now refer to a new, smaller vector. The chart of `subs`
2264: is `[min(subpointMap),max(subpointMap)+1)`
2266: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2268: Developer Notes:
2269: The use of the term Subdomain is unneeded and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`
2271: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2272: @*/
2273: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2274: {
2275: PetscFunctionBegin;
2276: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2277: PetscFunctionReturn(PETSC_SUCCESS);
2278: }
2280: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2281: {
2282: PetscInt p;
2283: PetscMPIInt rank;
2285: PetscFunctionBegin;
2286: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2287: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2288: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2289: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2290: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2291: PetscInt b;
2293: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2294: if (s->bcIndices) {
2295: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2296: }
2297: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2298: } else {
2299: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2300: }
2301: }
2302: PetscCall(PetscViewerFlush(viewer));
2303: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2304: if (s->sym) {
2305: PetscCall(PetscViewerASCIIPushTab(viewer));
2306: PetscCall(PetscSectionSymView(s->sym, viewer));
2307: PetscCall(PetscViewerASCIIPopTab(viewer));
2308: }
2309: PetscFunctionReturn(PETSC_SUCCESS);
2310: }
2312: /*@C
2313: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2315: Collective
2317: Input Parameters:
2318: + A - the `PetscSection` object to view
2319: . obj - Optional object that provides the options prefix used for the options
2320: - name - command line option
2322: Level: intermediate
2324: Note:
2325: See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`
2327: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2328: @*/
2329: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2330: {
2331: PetscFunctionBegin;
2333: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2334: PetscFunctionReturn(PETSC_SUCCESS);
2335: }
2337: /*@C
2338: PetscSectionView - Views a `PetscSection`
2340: Collective
2342: Input Parameters:
2343: + s - the `PetscSection` object to view
2344: - viewer - the viewer
2346: Level: beginner
2348: Note:
2349: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2350: distribution independent data, such as dofs, offsets, constraint dofs,
2351: and constraint indices. Points that have negative dofs, for instance,
2352: are not saved as they represent points owned by other processes.
2353: Point numbering and rank assignment is currently not stored.
2354: The saved section can be loaded with `PetscSectionLoad()`.
2356: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2357: @*/
2358: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2359: {
2360: PetscBool isascii, ishdf5;
2361: PetscInt f;
2363: PetscFunctionBegin;
2365: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2367: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2368: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2369: if (isascii) {
2370: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2371: if (s->numFields) {
2372: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2373: for (f = 0; f < s->numFields; ++f) {
2374: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2375: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2376: }
2377: } else {
2378: PetscCall(PetscSectionView_ASCII(s, viewer));
2379: }
2380: } else if (ishdf5) {
2381: #if PetscDefined(HAVE_HDF5)
2382: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2383: #else
2384: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2385: #endif
2386: }
2387: PetscFunctionReturn(PETSC_SUCCESS);
2388: }
2390: /*@C
2391: PetscSectionLoad - Loads a `PetscSection`
2393: Collective
2395: Input Parameters:
2396: + s - the `PetscSection` object to load
2397: - viewer - the viewer
2399: Level: beginner
2401: Note:
2402: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2403: a section saved with `PetscSectionView()`. The number of processes
2404: used here (N) does not need to be the same as that used when saving.
2405: After calling this function, the chart of s on rank i will be set
2406: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2407: saved section points.
2409: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2410: @*/
2411: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2412: {
2413: PetscBool ishdf5;
2415: PetscFunctionBegin;
2418: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2419: if (ishdf5) {
2420: #if PetscDefined(HAVE_HDF5)
2421: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2422: PetscFunctionReturn(PETSC_SUCCESS);
2423: #else
2424: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2425: #endif
2426: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2427: }
2429: /*@
2430: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2432: Input Parameter:
2433: . section - The `PetscSection`
2435: Level: intermediate
2437: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2438: @*/
2439: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2440: {
2441: PetscSectionClosurePermVal clVal;
2443: PetscFunctionBegin;
2444: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2445: kh_foreach_value(section->clHash, clVal, {
2446: PetscCall(PetscFree(clVal.perm));
2447: PetscCall(PetscFree(clVal.invPerm));
2448: });
2449: kh_destroy(ClPerm, section->clHash);
2450: section->clHash = NULL;
2451: PetscFunctionReturn(PETSC_SUCCESS);
2452: }
2454: /*@
2455: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2457: Not Collective
2459: Input Parameter:
2460: . s - the `PetscSection`
2462: Level: beginner
2464: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2465: @*/
2466: PetscErrorCode PetscSectionReset(PetscSection s)
2467: {
2468: PetscInt f, c;
2470: PetscFunctionBegin;
2472: for (f = 0; f < s->numFields; ++f) {
2473: PetscCall(PetscSectionDestroy(&s->field[f]));
2474: PetscCall(PetscFree(s->fieldNames[f]));
2475: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2476: PetscCall(PetscFree(s->compNames[f]));
2477: }
2478: PetscCall(PetscFree(s->numFieldComponents));
2479: PetscCall(PetscFree(s->fieldNames));
2480: PetscCall(PetscFree(s->compNames));
2481: PetscCall(PetscFree(s->field));
2482: PetscCall(PetscSectionDestroy(&s->bc));
2483: PetscCall(PetscFree(s->bcIndices));
2484: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2485: PetscCall(PetscSectionDestroy(&s->clSection));
2486: PetscCall(ISDestroy(&s->clPoints));
2487: PetscCall(ISDestroy(&s->perm));
2488: PetscCall(PetscSectionResetClosurePermutation(s));
2489: PetscCall(PetscSectionSymDestroy(&s->sym));
2490: PetscCall(PetscSectionDestroy(&s->clSection));
2491: PetscCall(ISDestroy(&s->clPoints));
2492: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2493: s->pStart = -1;
2494: s->pEnd = -1;
2495: s->maxDof = 0;
2496: s->setup = PETSC_FALSE;
2497: s->numFields = 0;
2498: s->clObj = NULL;
2499: PetscFunctionReturn(PETSC_SUCCESS);
2500: }
2502: /*@
2503: PetscSectionDestroy - Frees a `PetscSection`
2505: Not Collective
2507: Input Parameter:
2508: . s - the `PetscSection`
2510: Level: beginner
2512: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2513: @*/
2514: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2515: {
2516: PetscFunctionBegin;
2517: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2519: if (--((PetscObject)(*s))->refct > 0) {
2520: *s = NULL;
2521: PetscFunctionReturn(PETSC_SUCCESS);
2522: }
2523: PetscCall(PetscSectionReset(*s));
2524: PetscCall(PetscHeaderDestroy(s));
2525: PetscFunctionReturn(PETSC_SUCCESS);
2526: }
2528: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2529: {
2530: const PetscInt p = point - s->pStart;
2532: PetscFunctionBegin;
2534: *values = &baseArray[s->atlasOff[p]];
2535: PetscFunctionReturn(PETSC_SUCCESS);
2536: }
2538: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2539: {
2540: PetscInt *array;
2541: const PetscInt p = point - s->pStart;
2542: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2543: PetscInt cDim = 0;
2545: PetscFunctionBegin;
2547: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2548: array = &baseArray[s->atlasOff[p]];
2549: if (!cDim) {
2550: if (orientation >= 0) {
2551: const PetscInt dim = s->atlasDof[p];
2552: PetscInt i;
2554: if (mode == INSERT_VALUES) {
2555: for (i = 0; i < dim; ++i) array[i] = values[i];
2556: } else {
2557: for (i = 0; i < dim; ++i) array[i] += values[i];
2558: }
2559: } else {
2560: PetscInt offset = 0;
2561: PetscInt j = -1, field, i;
2563: for (field = 0; field < s->numFields; ++field) {
2564: const PetscInt dim = s->field[field]->atlasDof[p];
2566: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2567: offset += dim;
2568: }
2569: }
2570: } else {
2571: if (orientation >= 0) {
2572: const PetscInt dim = s->atlasDof[p];
2573: PetscInt cInd = 0, i;
2574: const PetscInt *cDof;
2576: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2577: if (mode == INSERT_VALUES) {
2578: for (i = 0; i < dim; ++i) {
2579: if ((cInd < cDim) && (i == cDof[cInd])) {
2580: ++cInd;
2581: continue;
2582: }
2583: array[i] = values[i];
2584: }
2585: } else {
2586: for (i = 0; i < dim; ++i) {
2587: if ((cInd < cDim) && (i == cDof[cInd])) {
2588: ++cInd;
2589: continue;
2590: }
2591: array[i] += values[i];
2592: }
2593: }
2594: } else {
2595: const PetscInt *cDof;
2596: PetscInt offset = 0;
2597: PetscInt cOffset = 0;
2598: PetscInt j = 0, field;
2600: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2601: for (field = 0; field < s->numFields; ++field) {
2602: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2603: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2604: const PetscInt sDim = dim - tDim;
2605: PetscInt cInd = 0, i, k;
2607: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2608: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2609: ++cInd;
2610: continue;
2611: }
2612: array[j] = values[k];
2613: }
2614: offset += dim;
2615: cOffset += dim - tDim;
2616: }
2617: }
2618: }
2619: PetscFunctionReturn(PETSC_SUCCESS);
2620: }
2622: /*@C
2623: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2625: Not Collective
2627: Input Parameter:
2628: . s - The `PetscSection`
2630: Output Parameter:
2631: . hasConstraints - flag indicating that the section has constrained dofs
2633: Level: intermediate
2635: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2636: @*/
2637: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2638: {
2639: PetscFunctionBegin;
2641: PetscAssertPointer(hasConstraints, 2);
2642: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2643: PetscFunctionReturn(PETSC_SUCCESS);
2644: }
2646: /*@C
2647: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2649: Not Collective
2651: Input Parameters:
2652: + s - The `PetscSection`
2653: - point - The point
2655: Output Parameter:
2656: . indices - The constrained dofs
2658: Level: intermediate
2660: Fortran Notes:
2661: Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2663: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2664: @*/
2665: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2666: {
2667: PetscFunctionBegin;
2669: if (s->bc) {
2670: PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2671: } else *indices = NULL;
2672: PetscFunctionReturn(PETSC_SUCCESS);
2673: }
2675: /*@C
2676: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2678: Not Collective
2680: Input Parameters:
2681: + s - The `PetscSection`
2682: . point - The point
2683: - indices - The constrained dofs
2685: Level: intermediate
2687: Fortran Notes:
2688: Use `PetscSectionSetConstraintIndicesF90()`
2690: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2691: @*/
2692: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2693: {
2694: PetscFunctionBegin;
2696: if (s->bc) {
2697: const PetscInt dof = s->atlasDof[point];
2698: const PetscInt cdof = s->bc->atlasDof[point];
2699: PetscInt d;
2701: for (d = 0; d < cdof; ++d) PetscCheck(indices[d] < dof, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " dof %" PetscInt_FMT ", invalid constraint index[%" PetscInt_FMT "]: %" PetscInt_FMT, point, dof, d, indices[d]);
2702: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2703: }
2704: PetscFunctionReturn(PETSC_SUCCESS);
2705: }
2707: /*@C
2708: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2710: Not Collective
2712: Input Parameters:
2713: + s - The `PetscSection`
2714: . field - The field number
2715: - point - The point
2717: Output Parameter:
2718: . indices - The constrained dofs sorted in ascending order
2720: Level: intermediate
2722: Note:
2723: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by `PetscSectionGetConstraintDof()`.
2725: Fortran Notes:
2726: Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2728: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2729: @*/
2730: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2731: {
2732: PetscFunctionBegin;
2734: PetscAssertPointer(indices, 4);
2735: PetscSectionCheckValidField(field, s->numFields);
2736: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2737: PetscFunctionReturn(PETSC_SUCCESS);
2738: }
2740: /*@C
2741: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2743: Not Collective
2745: Input Parameters:
2746: + s - The `PetscSection`
2747: . point - The point
2748: . field - The field number
2749: - indices - The constrained dofs
2751: Level: intermediate
2753: Fortran Notes:
2754: Use `PetscSectionSetFieldConstraintIndicesF90()`
2756: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2757: @*/
2758: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2759: {
2760: PetscFunctionBegin;
2762: if (PetscDefined(USE_DEBUG)) {
2763: PetscInt nfdof;
2765: PetscCall(PetscSectionGetFieldConstraintDof(s, point, field, &nfdof));
2766: if (nfdof) PetscAssertPointer(indices, 4);
2767: }
2768: PetscSectionCheckValidField(field, s->numFields);
2769: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2770: PetscFunctionReturn(PETSC_SUCCESS);
2771: }
2773: /*@
2774: PetscSectionPermute - Reorder the section according to the input point permutation
2776: Collective
2778: Input Parameters:
2779: + section - The `PetscSection` object
2780: - permutation - The point permutation, old point p becomes new point perm[p]
2782: Output Parameter:
2783: . sectionNew - The permuted `PetscSection`
2785: Level: intermediate
2787: Note:
2788: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
2790: Compare to `PetscSectionSetPermutation()`
2792: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2793: @*/
2794: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2795: {
2796: PetscSection s = section, sNew;
2797: const PetscInt *perm;
2798: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2800: PetscFunctionBegin;
2803: PetscAssertPointer(sectionNew, 3);
2804: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2805: PetscCall(PetscSectionGetNumFields(s, &numFields));
2806: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2807: for (f = 0; f < numFields; ++f) {
2808: const char *name;
2809: PetscInt numComp;
2811: PetscCall(PetscSectionGetFieldName(s, f, &name));
2812: PetscCall(PetscSectionSetFieldName(sNew, f, name));
2813: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2814: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2815: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2816: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2817: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2818: }
2819: }
2820: PetscCall(ISGetLocalSize(permutation, &numPoints));
2821: PetscCall(ISGetIndices(permutation, &perm));
2822: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2823: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2824: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2825: for (p = pStart; p < pEnd; ++p) {
2826: PetscInt dof, cdof;
2828: PetscCall(PetscSectionGetDof(s, p, &dof));
2829: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2830: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2831: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2832: for (f = 0; f < numFields; ++f) {
2833: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2834: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2835: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2836: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2837: }
2838: }
2839: PetscCall(PetscSectionSetUp(sNew));
2840: for (p = pStart; p < pEnd; ++p) {
2841: const PetscInt *cind;
2842: PetscInt cdof;
2844: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2845: if (cdof) {
2846: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
2847: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
2848: }
2849: for (f = 0; f < numFields; ++f) {
2850: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2851: if (cdof) {
2852: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
2853: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
2854: }
2855: }
2856: }
2857: PetscCall(ISRestoreIndices(permutation, &perm));
2858: *sectionNew = sNew;
2859: PetscFunctionReturn(PETSC_SUCCESS);
2860: }
2862: /*@
2863: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
2865: Collective
2867: Input Parameters:
2868: + section - The `PetscSection`
2869: . obj - A `PetscObject` which serves as the key for this index
2870: . clSection - `PetscSection` giving the size of the closure of each point
2871: - clPoints - `IS` giving the points in each closure
2873: Level: advanced
2875: Note:
2876: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
2878: Developer Notes:
2879: The information provided here is completely opaque
2881: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2882: @*/
2883: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2884: {
2885: PetscFunctionBegin;
2889: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
2890: section->clObj = obj;
2891: PetscCall(PetscObjectReference((PetscObject)clSection));
2892: PetscCall(PetscObjectReference((PetscObject)clPoints));
2893: PetscCall(PetscSectionDestroy(§ion->clSection));
2894: PetscCall(ISDestroy(§ion->clPoints));
2895: section->clSection = clSection;
2896: section->clPoints = clPoints;
2897: PetscFunctionReturn(PETSC_SUCCESS);
2898: }
2900: /*@
2901: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
2903: Collective
2905: Input Parameters:
2906: + section - The `PetscSection`
2907: - obj - A `PetscObject` which serves as the key for this index
2909: Output Parameters:
2910: + clSection - `PetscSection` giving the size of the closure of each point
2911: - clPoints - `IS` giving the points in each closure
2913: Level: advanced
2915: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2916: @*/
2917: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2918: {
2919: PetscFunctionBegin;
2920: if (section->clObj == obj) {
2921: if (clSection) *clSection = section->clSection;
2922: if (clPoints) *clPoints = section->clPoints;
2923: } else {
2924: if (clSection) *clSection = NULL;
2925: if (clPoints) *clPoints = NULL;
2926: }
2927: PetscFunctionReturn(PETSC_SUCCESS);
2928: }
2930: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2931: {
2932: PetscInt i;
2933: khiter_t iter;
2934: int new_entry;
2935: PetscSectionClosurePermKey key = {depth, clSize};
2936: PetscSectionClosurePermVal *val;
2938: PetscFunctionBegin;
2939: if (section->clObj != obj) {
2940: PetscCall(PetscSectionDestroy(§ion->clSection));
2941: PetscCall(ISDestroy(§ion->clPoints));
2942: }
2943: section->clObj = obj;
2944: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
2945: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2946: val = &kh_val(section->clHash, iter);
2947: if (!new_entry) {
2948: PetscCall(PetscFree(val->perm));
2949: PetscCall(PetscFree(val->invPerm));
2950: }
2951: if (mode == PETSC_COPY_VALUES) {
2952: PetscCall(PetscMalloc1(clSize, &val->perm));
2953: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
2954: } else if (mode == PETSC_OWN_POINTER) {
2955: val->perm = clPerm;
2956: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2957: PetscCall(PetscMalloc1(clSize, &val->invPerm));
2958: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2959: PetscFunctionReturn(PETSC_SUCCESS);
2960: }
2962: /*@
2963: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2965: Not Collective
2967: Input Parameters:
2968: + section - The `PetscSection`
2969: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2970: . depth - Depth of points on which to apply the given permutation
2971: - perm - Permutation of the cell dof closure
2973: Level: intermediate
2975: Notes:
2976: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2977: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2978: each topology and degree.
2980: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2981: exotic/enriched spaces on mixed topology meshes.
2983: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2984: @*/
2985: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2986: {
2987: const PetscInt *clPerm = NULL;
2988: PetscInt clSize = 0;
2990: PetscFunctionBegin;
2991: if (perm) {
2992: PetscCall(ISGetLocalSize(perm, &clSize));
2993: PetscCall(ISGetIndices(perm, &clPerm));
2994: }
2995: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
2996: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
2997: PetscFunctionReturn(PETSC_SUCCESS);
2998: }
3000: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3001: {
3002: PetscFunctionBegin;
3003: if (section->clObj == obj) {
3004: PetscSectionClosurePermKey k = {depth, size};
3005: PetscSectionClosurePermVal v;
3007: PetscCall(PetscClPermGet(section->clHash, k, &v));
3008: if (perm) *perm = v.perm;
3009: } else {
3010: if (perm) *perm = NULL;
3011: }
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: /*@
3016: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3018: Not Collective
3020: Input Parameters:
3021: + section - The `PetscSection`
3022: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3023: . depth - Depth stratum on which to obtain closure permutation
3024: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3026: Output Parameter:
3027: . perm - The dof closure permutation
3029: Level: intermediate
3031: Note:
3032: The user must destroy the `IS` that is returned.
3034: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3035: @*/
3036: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3037: {
3038: const PetscInt *clPerm = NULL;
3040: PetscFunctionBegin;
3041: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3042: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3043: PetscFunctionReturn(PETSC_SUCCESS);
3044: }
3046: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3047: {
3048: PetscFunctionBegin;
3049: if (section->clObj == obj && section->clHash) {
3050: PetscSectionClosurePermKey k = {depth, size};
3051: PetscSectionClosurePermVal v;
3052: PetscCall(PetscClPermGet(section->clHash, k, &v));
3053: if (perm) *perm = v.invPerm;
3054: } else {
3055: if (perm) *perm = NULL;
3056: }
3057: PetscFunctionReturn(PETSC_SUCCESS);
3058: }
3060: /*@
3061: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3063: Not Collective
3065: Input Parameters:
3066: + section - The `PetscSection`
3067: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3068: . depth - Depth stratum on which to obtain closure permutation
3069: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3071: Output Parameter:
3072: . perm - The dof closure permutation
3074: Level: intermediate
3076: Note:
3077: The user must destroy the `IS` that is returned.
3079: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3080: @*/
3081: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3082: {
3083: const PetscInt *clPerm = NULL;
3085: PetscFunctionBegin;
3086: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3087: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3088: PetscFunctionReturn(PETSC_SUCCESS);
3089: }
3091: /*@
3092: PetscSectionGetField - Get the `PetscSection` associated with a single field
3094: Input Parameters:
3095: + s - The `PetscSection`
3096: - field - The field number
3098: Output Parameter:
3099: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3101: Level: intermediate
3103: Note:
3104: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3106: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3107: @*/
3108: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3109: {
3110: PetscFunctionBegin;
3112: PetscAssertPointer(subs, 3);
3113: PetscSectionCheckValidField(field, s->numFields);
3114: *subs = s->field[field];
3115: PetscFunctionReturn(PETSC_SUCCESS);
3116: }
3118: PetscClassId PETSC_SECTION_SYM_CLASSID;
3119: PetscFunctionList PetscSectionSymList = NULL;
3121: /*@
3122: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3124: Collective
3126: Input Parameter:
3127: . comm - the MPI communicator
3129: Output Parameter:
3130: . sym - pointer to the new set of symmetries
3132: Level: developer
3134: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3135: @*/
3136: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3137: {
3138: PetscFunctionBegin;
3139: PetscAssertPointer(sym, 2);
3140: PetscCall(ISInitializePackage());
3141: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3142: PetscFunctionReturn(PETSC_SUCCESS);
3143: }
3145: /*@C
3146: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3148: Collective
3150: Input Parameters:
3151: + sym - The section symmetry object
3152: - method - The name of the section symmetry type
3154: Level: developer
3156: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3157: @*/
3158: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3159: {
3160: PetscErrorCode (*r)(PetscSectionSym);
3161: PetscBool match;
3163: PetscFunctionBegin;
3165: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3166: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3168: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3169: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3170: PetscTryTypeMethod(sym, destroy);
3171: sym->ops->destroy = NULL;
3173: PetscCall((*r)(sym));
3174: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3175: PetscFunctionReturn(PETSC_SUCCESS);
3176: }
3178: /*@C
3179: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3181: Not Collective
3183: Input Parameter:
3184: . sym - The section symmetry
3186: Output Parameter:
3187: . type - The index set type name
3189: Level: developer
3191: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3192: @*/
3193: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3194: {
3195: PetscFunctionBegin;
3197: PetscAssertPointer(type, 2);
3198: *type = ((PetscObject)sym)->type_name;
3199: PetscFunctionReturn(PETSC_SUCCESS);
3200: }
3202: /*@C
3203: PetscSectionSymRegister - Registers a new section symmetry implementation
3205: Not Collective
3207: Input Parameters:
3208: + sname - The name of a new user-defined creation routine
3209: - function - The creation routine itself
3211: Level: developer
3213: Notes:
3214: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3216: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3217: @*/
3218: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3219: {
3220: PetscFunctionBegin;
3221: PetscCall(ISInitializePackage());
3222: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3223: PetscFunctionReturn(PETSC_SUCCESS);
3224: }
3226: /*@
3227: PetscSectionSymDestroy - Destroys a section symmetry.
3229: Collective
3231: Input Parameter:
3232: . sym - the section symmetry
3234: Level: developer
3236: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3237: @*/
3238: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3239: {
3240: SymWorkLink link, next;
3242: PetscFunctionBegin;
3243: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3245: if (--((PetscObject)(*sym))->refct > 0) {
3246: *sym = NULL;
3247: PetscFunctionReturn(PETSC_SUCCESS);
3248: }
3249: if ((*sym)->ops->destroy) PetscCall((*(*sym)->ops->destroy)(*sym));
3250: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3251: for (link = (*sym)->workin; link; link = next) {
3252: PetscInt **perms = (PetscInt **)link->perms;
3253: PetscScalar **rots = (PetscScalar **)link->rots;
3254: PetscCall(PetscFree2(perms, rots));
3255: next = link->next;
3256: PetscCall(PetscFree(link));
3257: }
3258: (*sym)->workin = NULL;
3259: PetscCall(PetscHeaderDestroy(sym));
3260: PetscFunctionReturn(PETSC_SUCCESS);
3261: }
3263: /*@C
3264: PetscSectionSymView - Displays a section symmetry
3266: Collective
3268: Input Parameters:
3269: + sym - the index set
3270: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3272: Level: developer
3274: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3275: @*/
3276: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3277: {
3278: PetscFunctionBegin;
3280: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3282: PetscCheckSameComm(sym, 1, viewer, 2);
3283: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3284: PetscTryTypeMethod(sym, view, viewer);
3285: PetscFunctionReturn(PETSC_SUCCESS);
3286: }
3288: /*@
3289: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3291: Collective
3293: Input Parameters:
3294: + section - the section describing data layout
3295: - sym - the symmetry describing the affect of orientation on the access of the data
3297: Level: developer
3299: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3300: @*/
3301: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3302: {
3303: PetscFunctionBegin;
3305: PetscCall(PetscSectionSymDestroy(&(section->sym)));
3306: if (sym) {
3308: PetscCheckSameComm(section, 1, sym, 2);
3309: PetscCall(PetscObjectReference((PetscObject)sym));
3310: }
3311: section->sym = sym;
3312: PetscFunctionReturn(PETSC_SUCCESS);
3313: }
3315: /*@
3316: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3318: Not Collective
3320: Input Parameter:
3321: . section - the section describing data layout
3323: Output Parameter:
3324: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3326: Level: developer
3328: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3329: @*/
3330: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3331: {
3332: PetscFunctionBegin;
3334: *sym = section->sym;
3335: PetscFunctionReturn(PETSC_SUCCESS);
3336: }
3338: /*@
3339: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3341: Collective
3343: Input Parameters:
3344: + section - the section describing data layout
3345: . field - the field number
3346: - sym - the symmetry describing the affect of orientation on the access of the data
3348: Level: developer
3350: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3351: @*/
3352: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3353: {
3354: PetscFunctionBegin;
3356: PetscSectionCheckValidField(field, section->numFields);
3357: PetscCall(PetscSectionSetSym(section->field[field], sym));
3358: PetscFunctionReturn(PETSC_SUCCESS);
3359: }
3361: /*@
3362: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3364: Collective
3366: Input Parameters:
3367: + section - the section describing data layout
3368: - field - the field number
3370: Output Parameter:
3371: . sym - the symmetry describing the affect of orientation on the access of the data
3373: Level: developer
3375: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3376: @*/
3377: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3378: {
3379: PetscFunctionBegin;
3381: PetscSectionCheckValidField(field, section->numFields);
3382: *sym = section->field[field]->sym;
3383: PetscFunctionReturn(PETSC_SUCCESS);
3384: }
3386: /*@C
3387: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3389: Not Collective
3391: Input Parameters:
3392: + section - the section
3393: . numPoints - the number of points
3394: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3395: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3396: context, see `DMPlexGetConeOrientation()`).
3398: Output Parameters:
3399: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3400: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3401: identity).
3403: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3404: .vb
3405: const PetscInt **perms;
3406: const PetscScalar **rots;
3407: PetscInt lOffset;
3409: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3410: for (i = 0, lOffset = 0; i < numPoints; i++) {
3411: PetscInt point = points[2*i], dof, sOffset;
3412: const PetscInt *perm = perms ? perms[i] : NULL;
3413: const PetscScalar *rot = rots ? rots[i] : NULL;
3415: PetscSectionGetDof(section,point,&dof);
3416: PetscSectionGetOffset(section,point,&sOffset);
3418: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3419: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3420: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3421: lOffset += dof;
3422: }
3423: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3424: .ve
3426: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3427: .vb
3428: const PetscInt **perms;
3429: const PetscScalar **rots;
3430: PetscInt lOffset;
3432: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3433: for (i = 0, lOffset = 0; i < numPoints; i++) {
3434: PetscInt point = points[2*i], dof, sOffset;
3435: const PetscInt *perm = perms ? perms[i] : NULL;
3436: const PetscScalar *rot = rots ? rots[i] : NULL;
3438: PetscSectionGetDof(section,point,&dof);
3439: PetscSectionGetOffset(section,point,&sOff);
3441: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3442: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3443: offset += dof;
3444: }
3445: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3446: .ve
3448: Level: developer
3450: Notes:
3451: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3453: Use `PetscSectionRestorePointSyms()` when finished with the data
3455: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3456: @*/
3457: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3458: {
3459: PetscSectionSym sym;
3461: PetscFunctionBegin;
3463: if (numPoints) PetscAssertPointer(points, 3);
3464: if (perms) *perms = NULL;
3465: if (rots) *rots = NULL;
3466: sym = section->sym;
3467: if (sym && (perms || rots)) {
3468: SymWorkLink link;
3470: if (sym->workin) {
3471: link = sym->workin;
3472: sym->workin = sym->workin->next;
3473: } else {
3474: PetscCall(PetscNew(&link));
3475: }
3476: if (numPoints > link->numPoints) {
3477: PetscInt **perms = (PetscInt **)link->perms;
3478: PetscScalar **rots = (PetscScalar **)link->rots;
3479: PetscCall(PetscFree2(perms, rots));
3480: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3481: link->numPoints = numPoints;
3482: }
3483: link->next = sym->workout;
3484: sym->workout = link;
3485: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3486: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3487: PetscCall((*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots));
3488: if (perms) *perms = link->perms;
3489: if (rots) *rots = link->rots;
3490: }
3491: PetscFunctionReturn(PETSC_SUCCESS);
3492: }
3494: /*@C
3495: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3497: Not Collective
3499: Input Parameters:
3500: + section - the section
3501: . numPoints - the number of points
3502: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3503: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3504: context, see `DMPlexGetConeOrientation()`).
3505: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3506: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3508: Level: developer
3510: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3511: @*/
3512: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3513: {
3514: PetscSectionSym sym;
3516: PetscFunctionBegin;
3518: sym = section->sym;
3519: if (sym && (perms || rots)) {
3520: SymWorkLink *p, link;
3522: for (p = &sym->workout; (link = *p); p = &link->next) {
3523: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3524: *p = link->next;
3525: link->next = sym->workin;
3526: sym->workin = link;
3527: if (perms) *perms = NULL;
3528: if (rots) *rots = NULL;
3529: PetscFunctionReturn(PETSC_SUCCESS);
3530: }
3531: }
3532: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3533: }
3534: PetscFunctionReturn(PETSC_SUCCESS);
3535: }
3537: /*@C
3538: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3540: Not Collective
3542: Input Parameters:
3543: + section - the section
3544: . field - the field of the section
3545: . numPoints - the number of points
3546: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3547: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3548: context, see `DMPlexGetConeOrientation()`).
3550: Output Parameters:
3551: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3552: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3553: identity).
3555: Level: developer
3557: Notes:
3558: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3560: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3562: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3563: @*/
3564: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3565: {
3566: PetscFunctionBegin;
3568: PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3569: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3570: PetscFunctionReturn(PETSC_SUCCESS);
3571: }
3573: /*@C
3574: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3576: Not Collective
3578: Input Parameters:
3579: + section - the section
3580: . field - the field number
3581: . numPoints - the number of points
3582: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3583: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3584: context, see `DMPlexGetConeOrientation()`).
3585: . perms - The permutations for the given orientations: set to NULL at conclusion
3586: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3588: Level: developer
3590: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3591: @*/
3592: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3593: {
3594: PetscFunctionBegin;
3596: PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3597: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3598: PetscFunctionReturn(PETSC_SUCCESS);
3599: }
3601: /*@
3602: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3604: Not Collective
3606: Input Parameter:
3607: . sym - the `PetscSectionSym`
3609: Output Parameter:
3610: . nsym - the equivalent symmetries
3612: Level: developer
3614: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3615: @*/
3616: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3617: {
3618: PetscFunctionBegin;
3621: PetscTryTypeMethod(sym, copy, nsym);
3622: PetscFunctionReturn(PETSC_SUCCESS);
3623: }
3625: /*@
3626: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3628: Collective
3630: Input Parameters:
3631: + sym - the `PetscSectionSym`
3632: - migrationSF - the distribution map from roots to leaves
3634: Output Parameter:
3635: . dsym - the redistributed symmetries
3637: Level: developer
3639: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3640: @*/
3641: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3642: {
3643: PetscFunctionBegin;
3646: PetscAssertPointer(dsym, 3);
3647: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3648: PetscFunctionReturn(PETSC_SUCCESS);
3649: }
3651: /*@
3652: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3654: Not Collective
3656: Input Parameter:
3657: . s - the global `PetscSection`
3659: Output Parameter:
3660: . flg - the flag
3662: Level: developer
3664: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3665: @*/
3666: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3667: {
3668: PetscFunctionBegin;
3670: *flg = s->useFieldOff;
3671: PetscFunctionReturn(PETSC_SUCCESS);
3672: }
3674: /*@
3675: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3677: Not Collective
3679: Input Parameters:
3680: + s - the global `PetscSection`
3681: - flg - the flag
3683: Level: developer
3685: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3686: @*/
3687: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3688: {
3689: PetscFunctionBegin;
3691: s->useFieldOff = flg;
3692: PetscFunctionReturn(PETSC_SUCCESS);
3693: }
3695: #define PetscSectionExpandPoints_Loop(TYPE) \
3696: do { \
3697: PetscInt i, n, o0, o1, size; \
3698: TYPE *a0 = (TYPE *)origArray, *a1; \
3699: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3700: PetscCall(PetscMalloc1(size, &a1)); \
3701: for (i = 0; i < npoints; i++) { \
3702: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3703: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3704: PetscCall(PetscSectionGetDof(s, i, &n)); \
3705: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3706: } \
3707: *newArray = (void *)a1; \
3708: } while (0)
3710: /*@
3711: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3713: Not Collective
3715: Input Parameters:
3716: + origSection - the `PetscSection` describing the layout of the array
3717: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3718: . origArray - the array; its size must be equal to the storage size of `origSection`
3719: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
3721: Output Parameters:
3722: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3723: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
3725: Level: developer
3727: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3728: @*/
3729: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3730: {
3731: PetscSection s;
3732: const PetscInt *points_;
3733: PetscInt i, n, npoints, pStart, pEnd;
3734: PetscMPIInt unitsize;
3736: PetscFunctionBegin;
3738: PetscAssertPointer(origArray, 3);
3740: if (newSection) PetscAssertPointer(newSection, 5);
3741: if (newArray) PetscAssertPointer(newArray, 6);
3742: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3743: PetscCall(ISGetLocalSize(points, &npoints));
3744: PetscCall(ISGetIndices(points, &points_));
3745: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3746: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3747: PetscCall(PetscSectionSetChart(s, 0, npoints));
3748: for (i = 0; i < npoints; i++) {
3749: PetscCheck(points_[i] >= pStart && points_[i] < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " (index %" PetscInt_FMT ") in input IS out of input section's chart", points_[i], i);
3750: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3751: PetscCall(PetscSectionSetDof(s, i, n));
3752: }
3753: PetscCall(PetscSectionSetUp(s));
3754: if (newArray) {
3755: if (dataType == MPIU_INT) {
3756: PetscSectionExpandPoints_Loop(PetscInt);
3757: } else if (dataType == MPIU_SCALAR) {
3758: PetscSectionExpandPoints_Loop(PetscScalar);
3759: } else if (dataType == MPIU_REAL) {
3760: PetscSectionExpandPoints_Loop(PetscReal);
3761: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3762: }
3763: if (newSection) {
3764: *newSection = s;
3765: } else {
3766: PetscCall(PetscSectionDestroy(&s));
3767: }
3768: PetscCall(ISRestoreIndices(points, &points_));
3769: PetscFunctionReturn(PETSC_SUCCESS);
3770: }