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;
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 Note:
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 Note:
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;
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: . section - 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;
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;
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;
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;
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 Note:
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;
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 Note:
499: The function name should have Field in it since they are field components.
501: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
502: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
503: @*/
504: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
505: {
506: PetscFunctionBegin;
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 Note:
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;
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 Note:
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
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) {
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;
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 contiquously, 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 contiquously, 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;
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;
862: 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);
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
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: 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);
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;
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
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;
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;
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;
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;
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;
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;
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 Note:
1523: This is a terrible function name
1525: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
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;
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 Note:
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;
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, this 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;
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
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, this 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;
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
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;
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);
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;
1949: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
1950: PetscCall(PetscSectionSetFieldName(*subs, f, name));
1951: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
1952: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
1953: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1954: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
1955: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
1956: }
1957: }
1958: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1959: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
1960: for (p = pStart; p < pEnd; ++p) {
1961: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1963: for (f = 0; f < len; ++f) {
1964: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1965: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
1966: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1967: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
1968: dof += fdof;
1969: cdof += cfdof;
1970: }
1971: PetscCall(PetscSectionSetDof(*subs, p, dof));
1972: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
1973: maxCdof = PetscMax(cdof, maxCdof);
1974: }
1975: PetscCall(PetscSectionSetUp(*subs));
1976: if (maxCdof) {
1977: PetscInt *indices;
1979: PetscCall(PetscMalloc1(maxCdof, &indices));
1980: for (p = pStart; p < pEnd; ++p) {
1981: PetscInt cdof;
1983: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
1984: if (cdof) {
1985: const PetscInt *oldIndices = NULL;
1986: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1988: for (f = 0; f < len; ++f) {
1989: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1990: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1991: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
1992: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
1993: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1994: numConst += cfdof;
1995: fOff += fdof;
1996: }
1997: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
1998: }
1999: }
2000: PetscCall(PetscFree(indices));
2001: }
2002: PetscFunctionReturn(PETSC_SUCCESS);
2003: }
2005: /*@
2006: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2008: Collective
2010: Input Parameters:
2011: + s - the input sections
2012: - len - the number of input sections
2014: Output Parameter:
2015: . supers - the supersection
2017: Level: advanced
2019: Notes:
2020: The section offsets now refer to a new, larger vector.
2022: Developer Note:
2023: Needs to explain how the sections are composed
2025: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2026: @*/
2027: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2028: {
2029: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
2031: PetscFunctionBegin;
2032: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2033: for (i = 0; i < len; ++i) {
2034: PetscInt nf, pStarti, pEndi;
2036: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2037: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2038: pStart = PetscMin(pStart, pStarti);
2039: pEnd = PetscMax(pEnd, pEndi);
2040: Nf += nf;
2041: }
2042: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2043: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2044: for (i = 0, f = 0; i < len; ++i) {
2045: PetscInt nf, fi, ci;
2047: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2048: for (fi = 0; fi < nf; ++fi, ++f) {
2049: const char *name = NULL;
2050: PetscInt numComp = 0;
2052: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2053: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2054: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2055: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2056: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2057: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2058: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2059: }
2060: }
2061: }
2062: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2063: for (p = pStart; p < pEnd; ++p) {
2064: PetscInt dof = 0, cdof = 0;
2066: for (i = 0, f = 0; i < len; ++i) {
2067: PetscInt nf, fi, pStarti, pEndi;
2068: PetscInt fdof = 0, cfdof = 0;
2070: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2071: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2072: if ((p < pStarti) || (p >= pEndi)) continue;
2073: for (fi = 0; fi < nf; ++fi, ++f) {
2074: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2075: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2076: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2077: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2078: dof += fdof;
2079: cdof += cfdof;
2080: }
2081: }
2082: PetscCall(PetscSectionSetDof(*supers, p, dof));
2083: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2084: maxCdof = PetscMax(cdof, maxCdof);
2085: }
2086: PetscCall(PetscSectionSetUp(*supers));
2087: if (maxCdof) {
2088: PetscInt *indices;
2090: PetscCall(PetscMalloc1(maxCdof, &indices));
2091: for (p = pStart; p < pEnd; ++p) {
2092: PetscInt cdof;
2094: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2095: if (cdof) {
2096: PetscInt dof, numConst = 0, fOff = 0;
2098: for (i = 0, f = 0; i < len; ++i) {
2099: const PetscInt *oldIndices = NULL;
2100: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2102: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2103: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2104: if ((p < pStarti) || (p >= pEndi)) continue;
2105: for (fi = 0; fi < nf; ++fi, ++f) {
2106: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2107: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2108: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2109: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2110: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2111: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2112: numConst += cfdof;
2113: }
2114: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2115: fOff += dof;
2116: }
2117: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2118: }
2119: }
2120: PetscCall(PetscFree(indices));
2121: }
2122: PetscFunctionReturn(PETSC_SUCCESS);
2123: }
2125: PetscErrorCode PetscSectionCreateSubplexSection_Internal(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
2126: {
2127: const PetscInt *points = NULL, *indices = NULL;
2128: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2130: PetscFunctionBegin;
2134: PetscCall(PetscSectionGetNumFields(s, &numFields));
2135: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2136: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2137: for (f = 0; f < numFields; ++f) {
2138: const char *name = NULL;
2139: PetscInt numComp = 0;
2141: PetscCall(PetscSectionGetFieldName(s, f, &name));
2142: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2143: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2144: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2145: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2146: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2147: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2148: }
2149: }
2150: /* For right now, we do not try to squeeze the subchart */
2151: if (subpointMap) {
2152: PetscCall(ISGetSize(subpointMap, &numSubpoints));
2153: PetscCall(ISGetIndices(subpointMap, &points));
2154: }
2155: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2156: if (renumberPoints) {
2157: spStart = 0;
2158: spEnd = numSubpoints;
2159: } else {
2160: PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd));
2161: ++spEnd;
2162: }
2163: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2164: for (p = pStart; p < pEnd; ++p) {
2165: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2167: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2168: if (subp < 0) continue;
2169: if (!renumberPoints) subp = p;
2170: for (f = 0; f < numFields; ++f) {
2171: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2172: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2173: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2174: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2175: }
2176: PetscCall(PetscSectionGetDof(s, p, &dof));
2177: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2178: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2179: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2180: }
2181: PetscCall(PetscSectionSetUp(*subs));
2182: /* Change offsets to original offsets */
2183: for (p = pStart; p < pEnd; ++p) {
2184: PetscInt off, foff = 0;
2186: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2187: if (subp < 0) continue;
2188: if (!renumberPoints) subp = p;
2189: for (f = 0; f < numFields; ++f) {
2190: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2191: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2192: }
2193: PetscCall(PetscSectionGetOffset(s, p, &off));
2194: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2195: }
2196: /* Copy constraint indices */
2197: for (subp = spStart; subp < spEnd; ++subp) {
2198: PetscInt cdof;
2200: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2201: if (cdof) {
2202: for (f = 0; f < numFields; ++f) {
2203: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2204: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2205: }
2206: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2207: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2208: }
2209: }
2210: if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points));
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: /*@
2215: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2217: Collective
2219: Input Parameters:
2220: + s - the `PetscSection`
2221: - subpointMap - a sorted list of points in the original mesh which are in the submesh
2223: Output Parameter:
2224: . subs - the subsection
2226: Level: advanced
2228: Notes:
2229: 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))`
2231: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2233: Developer Note:
2234: 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`
2236: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2237: @*/
2238: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2239: {
2240: PetscFunctionBegin;
2241: PetscCall(PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs));
2242: PetscFunctionReturn(PETSC_SUCCESS);
2243: }
2245: /*@
2246: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2248: Collective
2250: Input Parameters:
2251: + s - the `PetscSection`
2252: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2254: Output Parameter:
2255: . subs - the subsection
2257: Level: advanced
2259: Notes:
2260: 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`
2261: is `[min(subpointMap),max(subpointMap)+1)`
2263: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2265: Developer Notes:
2266: 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`
2268: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2269: @*/
2270: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2271: {
2272: PetscFunctionBegin;
2273: PetscCall(PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs));
2274: PetscFunctionReturn(PETSC_SUCCESS);
2275: }
2277: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2278: {
2279: PetscInt p;
2280: PetscMPIInt rank;
2282: PetscFunctionBegin;
2283: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2284: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2285: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2286: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2287: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2288: PetscInt b;
2290: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2291: if (s->bcIndices) {
2292: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2293: }
2294: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2295: } else {
2296: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2297: }
2298: }
2299: PetscCall(PetscViewerFlush(viewer));
2300: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2301: if (s->sym) {
2302: PetscCall(PetscViewerASCIIPushTab(viewer));
2303: PetscCall(PetscSectionSymView(s->sym, viewer));
2304: PetscCall(PetscViewerASCIIPopTab(viewer));
2305: }
2306: PetscFunctionReturn(PETSC_SUCCESS);
2307: }
2309: /*@C
2310: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2312: Collective
2314: Input Parameters:
2315: + A - the `PetscSection` object to view
2316: . obj - Optional object that provides the options prefix used for the options
2317: - name - command line option
2319: Level: intermediate
2321: Note:
2322: See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`
2324: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2325: @*/
2326: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2327: {
2328: PetscFunctionBegin;
2330: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2331: PetscFunctionReturn(PETSC_SUCCESS);
2332: }
2334: /*@C
2335: PetscSectionView - Views a `PetscSection`
2337: Collective
2339: Input Parameters:
2340: + s - the `PetscSection` object to view
2341: - v - the viewer
2343: Level: beginner
2345: Note:
2346: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2347: distribution independent data, such as dofs, offsets, constraint dofs,
2348: and constraint indices. Points that have negative dofs, for instance,
2349: are not saved as they represent points owned by other processes.
2350: Point numbering and rank assignment is currently not stored.
2351: The saved section can be loaded with `PetscSectionLoad()`.
2353: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2354: @*/
2355: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2356: {
2357: PetscBool isascii, ishdf5;
2358: PetscInt f;
2360: PetscFunctionBegin;
2362: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2364: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2365: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2366: if (isascii) {
2367: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2368: if (s->numFields) {
2369: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2370: for (f = 0; f < s->numFields; ++f) {
2371: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2372: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2373: }
2374: } else {
2375: PetscCall(PetscSectionView_ASCII(s, viewer));
2376: }
2377: } else if (ishdf5) {
2378: #if PetscDefined(HAVE_HDF5)
2379: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2380: #else
2381: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2382: #endif
2383: }
2384: PetscFunctionReturn(PETSC_SUCCESS);
2385: }
2387: /*@C
2388: PetscSectionLoad - Loads a `PetscSection`
2390: Collective
2392: Input Parameters:
2393: + s - the `PetscSection` object to load
2394: - v - the viewer
2396: Level: beginner
2398: Note:
2399: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2400: a section saved with `PetscSectionView()`. The number of processes
2401: used here (N) does not need to be the same as that used when saving.
2402: After calling this function, the chart of s on rank i will be set
2403: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2404: saved section points.
2406: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2407: @*/
2408: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2409: {
2410: PetscBool ishdf5;
2412: PetscFunctionBegin;
2415: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2416: if (ishdf5) {
2417: #if PetscDefined(HAVE_HDF5)
2418: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2419: PetscFunctionReturn(PETSC_SUCCESS);
2420: #else
2421: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2422: #endif
2423: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2424: }
2426: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2427: {
2428: PetscSectionClosurePermVal clVal;
2430: PetscFunctionBegin;
2431: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2432: kh_foreach_value(section->clHash, clVal, {
2433: PetscCall(PetscFree(clVal.perm));
2434: PetscCall(PetscFree(clVal.invPerm));
2435: });
2436: kh_destroy(ClPerm, section->clHash);
2437: section->clHash = NULL;
2438: PetscFunctionReturn(PETSC_SUCCESS);
2439: }
2441: /*@
2442: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2444: Not Collective
2446: Input Parameter:
2447: . s - the `PetscSection`
2449: Level: beginner
2451: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2452: @*/
2453: PetscErrorCode PetscSectionReset(PetscSection s)
2454: {
2455: PetscInt f, c;
2457: PetscFunctionBegin;
2459: for (f = 0; f < s->numFields; ++f) {
2460: PetscCall(PetscSectionDestroy(&s->field[f]));
2461: PetscCall(PetscFree(s->fieldNames[f]));
2462: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2463: PetscCall(PetscFree(s->compNames[f]));
2464: }
2465: PetscCall(PetscFree(s->numFieldComponents));
2466: PetscCall(PetscFree(s->fieldNames));
2467: PetscCall(PetscFree(s->compNames));
2468: PetscCall(PetscFree(s->field));
2469: PetscCall(PetscSectionDestroy(&s->bc));
2470: PetscCall(PetscFree(s->bcIndices));
2471: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2472: PetscCall(PetscSectionDestroy(&s->clSection));
2473: PetscCall(ISDestroy(&s->clPoints));
2474: PetscCall(ISDestroy(&s->perm));
2475: PetscCall(PetscSectionResetClosurePermutation(s));
2476: PetscCall(PetscSectionSymDestroy(&s->sym));
2477: PetscCall(PetscSectionDestroy(&s->clSection));
2478: PetscCall(ISDestroy(&s->clPoints));
2479: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2480: s->pStart = -1;
2481: s->pEnd = -1;
2482: s->maxDof = 0;
2483: s->setup = PETSC_FALSE;
2484: s->numFields = 0;
2485: s->clObj = NULL;
2486: PetscFunctionReturn(PETSC_SUCCESS);
2487: }
2489: /*@
2490: PetscSectionDestroy - Frees a `PetscSection`
2492: Not Collective
2494: Input Parameter:
2495: . s - the `PetscSection`
2497: Level: beginner
2499: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2500: @*/
2501: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2502: {
2503: PetscFunctionBegin;
2504: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2506: if (--((PetscObject)(*s))->refct > 0) {
2507: *s = NULL;
2508: PetscFunctionReturn(PETSC_SUCCESS);
2509: }
2510: PetscCall(PetscSectionReset(*s));
2511: PetscCall(PetscHeaderDestroy(s));
2512: PetscFunctionReturn(PETSC_SUCCESS);
2513: }
2515: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2516: {
2517: const PetscInt p = point - s->pStart;
2519: PetscFunctionBegin;
2521: *values = &baseArray[s->atlasOff[p]];
2522: PetscFunctionReturn(PETSC_SUCCESS);
2523: }
2525: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2526: {
2527: PetscInt *array;
2528: const PetscInt p = point - s->pStart;
2529: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2530: PetscInt cDim = 0;
2532: PetscFunctionBegin;
2534: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2535: array = &baseArray[s->atlasOff[p]];
2536: if (!cDim) {
2537: if (orientation >= 0) {
2538: const PetscInt dim = s->atlasDof[p];
2539: PetscInt i;
2541: if (mode == INSERT_VALUES) {
2542: for (i = 0; i < dim; ++i) array[i] = values[i];
2543: } else {
2544: for (i = 0; i < dim; ++i) array[i] += values[i];
2545: }
2546: } else {
2547: PetscInt offset = 0;
2548: PetscInt j = -1, field, i;
2550: for (field = 0; field < s->numFields; ++field) {
2551: const PetscInt dim = s->field[field]->atlasDof[p];
2553: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2554: offset += dim;
2555: }
2556: }
2557: } else {
2558: if (orientation >= 0) {
2559: const PetscInt dim = s->atlasDof[p];
2560: PetscInt cInd = 0, i;
2561: const PetscInt *cDof;
2563: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2564: if (mode == INSERT_VALUES) {
2565: for (i = 0; i < dim; ++i) {
2566: if ((cInd < cDim) && (i == cDof[cInd])) {
2567: ++cInd;
2568: continue;
2569: }
2570: array[i] = values[i];
2571: }
2572: } else {
2573: for (i = 0; i < dim; ++i) {
2574: if ((cInd < cDim) && (i == cDof[cInd])) {
2575: ++cInd;
2576: continue;
2577: }
2578: array[i] += values[i];
2579: }
2580: }
2581: } else {
2582: const PetscInt *cDof;
2583: PetscInt offset = 0;
2584: PetscInt cOffset = 0;
2585: PetscInt j = 0, field;
2587: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2588: for (field = 0; field < s->numFields; ++field) {
2589: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2590: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2591: const PetscInt sDim = dim - tDim;
2592: PetscInt cInd = 0, i, k;
2594: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2595: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2596: ++cInd;
2597: continue;
2598: }
2599: array[j] = values[k];
2600: }
2601: offset += dim;
2602: cOffset += dim - tDim;
2603: }
2604: }
2605: }
2606: PetscFunctionReturn(PETSC_SUCCESS);
2607: }
2609: /*@C
2610: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2612: Not Collective
2614: Input Parameter:
2615: . s - The `PetscSection`
2617: Output Parameter:
2618: . hasConstraints - flag indicating that the section has constrained dofs
2620: Level: intermediate
2622: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2623: @*/
2624: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2625: {
2626: PetscFunctionBegin;
2629: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2630: PetscFunctionReturn(PETSC_SUCCESS);
2631: }
2633: /*@C
2634: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2636: Not Collective
2638: Input Parameters:
2639: + s - The `PetscSection`
2640: - point - The point
2642: Output Parameter:
2643: . indices - The constrained dofs
2645: Level: intermediate
2647: Fortran Note:
2648: Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2650: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2651: @*/
2652: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2653: {
2654: PetscFunctionBegin;
2656: if (s->bc) {
2657: PetscCall(VecIntGetValuesSection(s->bcIndices, s->bc, point, indices));
2658: } else *indices = NULL;
2659: PetscFunctionReturn(PETSC_SUCCESS);
2660: }
2662: /*@C
2663: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2665: Not Collective
2667: Input Parameters:
2668: + s - The `PetscSection`
2669: . point - The point
2670: - indices - The constrained dofs
2672: Level: intermediate
2674: Fortran Note:
2675: Use `PetscSectionSetConstraintIndicesF90()`
2677: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2678: @*/
2679: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2680: {
2681: PetscFunctionBegin;
2683: if (s->bc) {
2684: const PetscInt dof = s->atlasDof[point];
2685: const PetscInt cdof = s->bc->atlasDof[point];
2686: PetscInt d;
2688: 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]);
2689: PetscCall(VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2690: }
2691: PetscFunctionReturn(PETSC_SUCCESS);
2692: }
2694: /*@C
2695: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2697: Not Collective
2699: Input Parameters:
2700: + s - The `PetscSection`
2701: . field - The field number
2702: - point - The point
2704: Output Parameter:
2705: . indices - The constrained dofs sorted in ascending order
2707: Level: intermediate
2709: Note:
2710: 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()`.
2712: Fortran Note:
2713: Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2715: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2716: @*/
2717: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2718: {
2719: PetscFunctionBegin;
2722: PetscSectionCheckValidField(field, s->numFields);
2723: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2724: PetscFunctionReturn(PETSC_SUCCESS);
2725: }
2727: /*@C
2728: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2730: Not Collective
2732: Input Parameters:
2733: + s - The `PetscSection`
2734: . point - The point
2735: . field - The field number
2736: - indices - The constrained dofs
2738: Level: intermediate
2740: Fortran Note:
2741: Use `PetscSectionSetFieldConstraintIndicesF90()`
2743: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2744: @*/
2745: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2746: {
2747: PetscFunctionBegin;
2749: if (PetscDefined(USE_DEBUG)) {
2750: PetscInt nfdof;
2752: PetscCall(PetscSectionGetFieldConstraintDof(s, point, field, &nfdof));
2754: }
2755: PetscSectionCheckValidField(field, s->numFields);
2756: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2757: PetscFunctionReturn(PETSC_SUCCESS);
2758: }
2760: /*@
2761: PetscSectionPermute - Reorder the section according to the input point permutation
2763: Collective
2765: Input Parameters:
2766: + section - The `PetscSection` object
2767: - perm - The point permutation, old point p becomes new point perm[p]
2769: Output Parameter:
2770: . sectionNew - The permuted `PetscSection`
2772: Level: intermediate
2774: Note:
2775: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
2777: Compare to `PetscSectionSetPermutation()`
2779: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2780: @*/
2781: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2782: {
2783: PetscSection s = section, sNew;
2784: const PetscInt *perm;
2785: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2787: PetscFunctionBegin;
2791: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2792: PetscCall(PetscSectionGetNumFields(s, &numFields));
2793: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2794: for (f = 0; f < numFields; ++f) {
2795: const char *name;
2796: PetscInt numComp;
2798: PetscCall(PetscSectionGetFieldName(s, f, &name));
2799: PetscCall(PetscSectionSetFieldName(sNew, f, name));
2800: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2801: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2802: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2803: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2804: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2805: }
2806: }
2807: PetscCall(ISGetLocalSize(permutation, &numPoints));
2808: PetscCall(ISGetIndices(permutation, &perm));
2809: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2810: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2811: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2812: for (p = pStart; p < pEnd; ++p) {
2813: PetscInt dof, cdof;
2815: PetscCall(PetscSectionGetDof(s, p, &dof));
2816: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2817: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2818: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2819: for (f = 0; f < numFields; ++f) {
2820: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2821: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2822: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2823: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2824: }
2825: }
2826: PetscCall(PetscSectionSetUp(sNew));
2827: for (p = pStart; p < pEnd; ++p) {
2828: const PetscInt *cind;
2829: PetscInt cdof;
2831: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2832: if (cdof) {
2833: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
2834: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
2835: }
2836: for (f = 0; f < numFields; ++f) {
2837: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2838: if (cdof) {
2839: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
2840: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
2841: }
2842: }
2843: }
2844: PetscCall(ISRestoreIndices(permutation, &perm));
2845: *sectionNew = sNew;
2846: PetscFunctionReturn(PETSC_SUCCESS);
2847: }
2849: /*@
2850: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
2852: Collective
2854: Input Parameters:
2855: + section - The `PetscSection`
2856: . obj - A `PetscObject` which serves as the key for this index
2857: . clSection - `PetscSection` giving the size of the closure of each point
2858: - clPoints - `IS` giving the points in each closure
2860: Level: advanced
2862: Note:
2863: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
2865: Developer Note:
2866: The information provided here is completely opaque
2868: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2869: @*/
2870: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2871: {
2872: PetscFunctionBegin;
2876: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
2877: section->clObj = obj;
2878: PetscCall(PetscObjectReference((PetscObject)clSection));
2879: PetscCall(PetscObjectReference((PetscObject)clPoints));
2880: PetscCall(PetscSectionDestroy(§ion->clSection));
2881: PetscCall(ISDestroy(§ion->clPoints));
2882: section->clSection = clSection;
2883: section->clPoints = clPoints;
2884: PetscFunctionReturn(PETSC_SUCCESS);
2885: }
2887: /*@
2888: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
2890: Collective
2892: Input Parameters:
2893: + section - The `PetscSection`
2894: - obj - A `PetscObject` which serves as the key for this index
2896: Output Parameters:
2897: + clSection - `PetscSection` giving the size of the closure of each point
2898: - clPoints - `IS` giving the points in each closure
2900: Level: advanced
2902: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2903: @*/
2904: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2905: {
2906: PetscFunctionBegin;
2907: if (section->clObj == obj) {
2908: if (clSection) *clSection = section->clSection;
2909: if (clPoints) *clPoints = section->clPoints;
2910: } else {
2911: if (clSection) *clSection = NULL;
2912: if (clPoints) *clPoints = NULL;
2913: }
2914: PetscFunctionReturn(PETSC_SUCCESS);
2915: }
2917: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2918: {
2919: PetscInt i;
2920: khiter_t iter;
2921: int new_entry;
2922: PetscSectionClosurePermKey key = {depth, clSize};
2923: PetscSectionClosurePermVal *val;
2925: PetscFunctionBegin;
2926: if (section->clObj != obj) {
2927: PetscCall(PetscSectionDestroy(§ion->clSection));
2928: PetscCall(ISDestroy(§ion->clPoints));
2929: }
2930: section->clObj = obj;
2931: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
2932: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2933: val = &kh_val(section->clHash, iter);
2934: if (!new_entry) {
2935: PetscCall(PetscFree(val->perm));
2936: PetscCall(PetscFree(val->invPerm));
2937: }
2938: if (mode == PETSC_COPY_VALUES) {
2939: PetscCall(PetscMalloc1(clSize, &val->perm));
2940: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
2941: } else if (mode == PETSC_OWN_POINTER) {
2942: val->perm = clPerm;
2943: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2944: PetscCall(PetscMalloc1(clSize, &val->invPerm));
2945: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2946: PetscFunctionReturn(PETSC_SUCCESS);
2947: }
2949: /*@
2950: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2952: Not Collective
2954: Input Parameters:
2955: + section - The `PetscSection`
2956: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2957: . depth - Depth of points on which to apply the given permutation
2958: - perm - Permutation of the cell dof closure
2960: Level: intermediate
2962: Notes:
2963: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2964: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2965: each topology and degree.
2967: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2968: exotic/enriched spaces on mixed topology meshes.
2970: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2971: @*/
2972: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2973: {
2974: const PetscInt *clPerm = NULL;
2975: PetscInt clSize = 0;
2977: PetscFunctionBegin;
2978: if (perm) {
2979: PetscCall(ISGetLocalSize(perm, &clSize));
2980: PetscCall(ISGetIndices(perm, &clPerm));
2981: }
2982: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
2983: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
2984: PetscFunctionReturn(PETSC_SUCCESS);
2985: }
2987: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2988: {
2989: PetscFunctionBegin;
2990: if (section->clObj == obj) {
2991: PetscSectionClosurePermKey k = {depth, size};
2992: PetscSectionClosurePermVal v;
2993: PetscCall(PetscClPermGet(section->clHash, k, &v));
2994: if (perm) *perm = v.perm;
2995: } else {
2996: if (perm) *perm = NULL;
2997: }
2998: PetscFunctionReturn(PETSC_SUCCESS);
2999: }
3001: /*@
3002: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3004: Not Collective
3006: Input Parameters:
3007: + section - The `PetscSection`
3008: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3009: . depth - Depth stratum on which to obtain closure permutation
3010: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3012: Output Parameter:
3013: . perm - The dof closure permutation
3015: Level: intermediate
3017: Note:
3018: The user must destroy the `IS` that is returned.
3020: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3021: @*/
3022: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3023: {
3024: const PetscInt *clPerm;
3026: PetscFunctionBegin;
3027: PetscCall(PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm));
3028: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3029: PetscFunctionReturn(PETSC_SUCCESS);
3030: }
3032: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3033: {
3034: PetscFunctionBegin;
3035: if (section->clObj == obj && section->clHash) {
3036: PetscSectionClosurePermKey k = {depth, size};
3037: PetscSectionClosurePermVal v;
3038: PetscCall(PetscClPermGet(section->clHash, k, &v));
3039: if (perm) *perm = v.invPerm;
3040: } else {
3041: if (perm) *perm = NULL;
3042: }
3043: PetscFunctionReturn(PETSC_SUCCESS);
3044: }
3046: /*@
3047: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3049: Not Collective
3051: Input Parameters:
3052: + section - The `PetscSection`
3053: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3054: . depth - Depth stratum on which to obtain closure permutation
3055: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3057: Output Parameter:
3058: . perm - The dof closure permutation
3060: Level: intermediate
3062: Note:
3063: The user must destroy the `IS` that is returned.
3065: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3066: @*/
3067: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3068: {
3069: const PetscInt *clPerm;
3071: PetscFunctionBegin;
3072: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3073: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3074: PetscFunctionReturn(PETSC_SUCCESS);
3075: }
3077: /*@
3078: PetscSectionGetField - Get the `PetscSection` associated with a single field
3080: Input Parameters:
3081: + s - The `PetscSection`
3082: - field - The field number
3084: Output Parameter:
3085: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3087: Level: intermediate
3089: Note:
3090: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3092: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3093: @*/
3094: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3095: {
3096: PetscFunctionBegin;
3099: PetscSectionCheckValidField(field, s->numFields);
3100: *subs = s->field[field];
3101: PetscFunctionReturn(PETSC_SUCCESS);
3102: }
3104: PetscClassId PETSC_SECTION_SYM_CLASSID;
3105: PetscFunctionList PetscSectionSymList = NULL;
3107: /*@
3108: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3110: Collective
3112: Input Parameter:
3113: . comm - the MPI communicator
3115: Output Parameter:
3116: . sym - pointer to the new set of symmetries
3118: Level: developer
3120: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3121: @*/
3122: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3123: {
3124: PetscFunctionBegin;
3126: PetscCall(ISInitializePackage());
3127: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3128: PetscFunctionReturn(PETSC_SUCCESS);
3129: }
3131: /*@C
3132: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3134: Collective
3136: Input Parameters:
3137: + sym - The section symmetry object
3138: - method - The name of the section symmetry type
3140: Level: developer
3142: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3143: @*/
3144: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3145: {
3146: PetscErrorCode (*r)(PetscSectionSym);
3147: PetscBool match;
3149: PetscFunctionBegin;
3151: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3152: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3154: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3155: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3156: PetscTryTypeMethod(sym, destroy);
3157: sym->ops->destroy = NULL;
3159: PetscCall((*r)(sym));
3160: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3161: PetscFunctionReturn(PETSC_SUCCESS);
3162: }
3164: /*@C
3165: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3167: Not Collective
3169: Input Parameter:
3170: . sym - The section symmetry
3172: Output Parameter:
3173: . type - The index set type name
3175: Level: developer
3177: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3178: @*/
3179: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3180: {
3181: PetscFunctionBegin;
3184: *type = ((PetscObject)sym)->type_name;
3185: PetscFunctionReturn(PETSC_SUCCESS);
3186: }
3188: /*@C
3189: PetscSectionSymRegister - Registers a new section symmetry implementation
3191: Not Collective
3193: Input Parameters:
3194: + sname - The name of a new user-defined creation routine
3195: - function - The creation routine itself
3197: Level: developer
3199: Notes:
3200: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3202: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3203: @*/
3204: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3205: {
3206: PetscFunctionBegin;
3207: PetscCall(ISInitializePackage());
3208: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3209: PetscFunctionReturn(PETSC_SUCCESS);
3210: }
3212: /*@
3213: PetscSectionSymDestroy - Destroys a section symmetry.
3215: Collective
3217: Input Parameter:
3218: . sym - the section symmetry
3220: Level: developer
3222: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
3223: @*/
3224: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3225: {
3226: SymWorkLink link, next;
3228: PetscFunctionBegin;
3229: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3231: if (--((PetscObject)(*sym))->refct > 0) {
3232: *sym = NULL;
3233: PetscFunctionReturn(PETSC_SUCCESS);
3234: }
3235: if ((*sym)->ops->destroy) PetscCall((*(*sym)->ops->destroy)(*sym));
3236: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3237: for (link = (*sym)->workin; link; link = next) {
3238: PetscInt **perms = (PetscInt **)link->perms;
3239: PetscScalar **rots = (PetscScalar **)link->rots;
3240: PetscCall(PetscFree2(perms, rots));
3241: next = link->next;
3242: PetscCall(PetscFree(link));
3243: }
3244: (*sym)->workin = NULL;
3245: PetscCall(PetscHeaderDestroy(sym));
3246: PetscFunctionReturn(PETSC_SUCCESS);
3247: }
3249: /*@C
3250: PetscSectionSymView - Displays a section symmetry
3252: Collective
3254: Input Parameters:
3255: + sym - the index set
3256: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3258: Level: developer
3260: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3261: @*/
3262: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3263: {
3264: PetscFunctionBegin;
3266: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3268: PetscCheckSameComm(sym, 1, viewer, 2);
3269: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3270: PetscTryTypeMethod(sym, view, viewer);
3271: PetscFunctionReturn(PETSC_SUCCESS);
3272: }
3274: /*@
3275: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3277: Collective
3279: Input Parameters:
3280: + section - the section describing data layout
3281: - sym - the symmetry describing the affect of orientation on the access of the data
3283: Level: developer
3285: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3286: @*/
3287: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3288: {
3289: PetscFunctionBegin;
3291: PetscCall(PetscSectionSymDestroy(&(section->sym)));
3292: if (sym) {
3294: PetscCheckSameComm(section, 1, sym, 2);
3295: PetscCall(PetscObjectReference((PetscObject)sym));
3296: }
3297: section->sym = sym;
3298: PetscFunctionReturn(PETSC_SUCCESS);
3299: }
3301: /*@
3302: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3304: Not Collective
3306: Input Parameter:
3307: . section - the section describing data layout
3309: Output Parameter:
3310: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3312: Level: developer
3314: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3315: @*/
3316: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3317: {
3318: PetscFunctionBegin;
3320: *sym = section->sym;
3321: PetscFunctionReturn(PETSC_SUCCESS);
3322: }
3324: /*@
3325: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3327: Collective
3329: Input Parameters:
3330: + section - the section describing data layout
3331: . field - the field number
3332: - sym - the symmetry describing the affect of orientation on the access of the data
3334: Level: developer
3336: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3337: @*/
3338: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3339: {
3340: PetscFunctionBegin;
3342: PetscSectionCheckValidField(field, section->numFields);
3343: PetscCall(PetscSectionSetSym(section->field[field], sym));
3344: PetscFunctionReturn(PETSC_SUCCESS);
3345: }
3347: /*@
3348: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3350: Collective
3352: Input Parameters:
3353: + section - the section describing data layout
3354: - field - the field number
3356: Output Parameter:
3357: . sym - the symmetry describing the affect of orientation on the access of the data
3359: Level: developer
3361: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3362: @*/
3363: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3364: {
3365: PetscFunctionBegin;
3367: PetscSectionCheckValidField(field, section->numFields);
3368: *sym = section->field[field]->sym;
3369: PetscFunctionReturn(PETSC_SUCCESS);
3370: }
3372: /*@C
3373: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3375: Not Collective
3377: Input Parameters:
3378: + section - the section
3379: . numPoints - the number of points
3380: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3381: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3382: context, see `DMPlexGetConeOrientation()`).
3384: Output Parameters:
3385: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3386: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3387: identity).
3389: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3390: .vb
3391: const PetscInt **perms;
3392: const PetscScalar **rots;
3393: PetscInt lOffset;
3395: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3396: for (i = 0, lOffset = 0; i < numPoints; i++) {
3397: PetscInt point = points[2*i], dof, sOffset;
3398: const PetscInt *perm = perms ? perms[i] : NULL;
3399: const PetscScalar *rot = rots ? rots[i] : NULL;
3401: PetscSectionGetDof(section,point,&dof);
3402: PetscSectionGetOffset(section,point,&sOffset);
3404: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3405: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3406: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3407: lOffset += dof;
3408: }
3409: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3410: .ve
3412: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3413: .vb
3414: const PetscInt **perms;
3415: const PetscScalar **rots;
3416: PetscInt lOffset;
3418: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3419: for (i = 0, lOffset = 0; i < numPoints; i++) {
3420: PetscInt point = points[2*i], dof, sOffset;
3421: const PetscInt *perm = perms ? perms[i] : NULL;
3422: const PetscScalar *rot = rots ? rots[i] : NULL;
3424: PetscSectionGetDof(section,point,&dof);
3425: PetscSectionGetOffset(section,point,&sOff);
3427: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3428: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3429: offset += dof;
3430: }
3431: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3432: .ve
3434: Level: developer
3436: Notes:
3437: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3439: Use `PetscSectionRestorePointSyms()` when finished with the data
3441: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3442: @*/
3443: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3444: {
3445: PetscSectionSym sym;
3447: PetscFunctionBegin;
3450: if (perms) *perms = NULL;
3451: if (rots) *rots = NULL;
3452: sym = section->sym;
3453: if (sym && (perms || rots)) {
3454: SymWorkLink link;
3456: if (sym->workin) {
3457: link = sym->workin;
3458: sym->workin = sym->workin->next;
3459: } else {
3460: PetscCall(PetscNew(&link));
3461: }
3462: if (numPoints > link->numPoints) {
3463: PetscInt **perms = (PetscInt **)link->perms;
3464: PetscScalar **rots = (PetscScalar **)link->rots;
3465: PetscCall(PetscFree2(perms, rots));
3466: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3467: link->numPoints = numPoints;
3468: }
3469: link->next = sym->workout;
3470: sym->workout = link;
3471: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3472: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3473: PetscCall((*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots));
3474: if (perms) *perms = link->perms;
3475: if (rots) *rots = link->rots;
3476: }
3477: PetscFunctionReturn(PETSC_SUCCESS);
3478: }
3480: /*@C
3481: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3483: Not Collective
3485: Input Parameters:
3486: + section - the section
3487: . numPoints - the number of points
3488: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3489: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3490: context, see `DMPlexGetConeOrientation()`).
3491: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3492: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3494: Level: developer
3496: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3497: @*/
3498: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3499: {
3500: PetscSectionSym sym;
3502: PetscFunctionBegin;
3504: sym = section->sym;
3505: if (sym && (perms || rots)) {
3506: SymWorkLink *p, link;
3508: for (p = &sym->workout; (link = *p); p = &link->next) {
3509: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3510: *p = link->next;
3511: link->next = sym->workin;
3512: sym->workin = link;
3513: if (perms) *perms = NULL;
3514: if (rots) *rots = NULL;
3515: PetscFunctionReturn(PETSC_SUCCESS);
3516: }
3517: }
3518: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3519: }
3520: PetscFunctionReturn(PETSC_SUCCESS);
3521: }
3523: /*@C
3524: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3526: Not Collective
3528: Input Parameters:
3529: + section - the section
3530: . field - the field of the section
3531: . numPoints - the number of points
3532: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3533: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3534: context, see `DMPlexGetConeOrientation()`).
3536: Output Parameters:
3537: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3538: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3539: identity).
3541: Level: developer
3543: Notes:
3544: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3546: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3548: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3549: @*/
3550: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3551: {
3552: PetscFunctionBegin;
3554: 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);
3555: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3556: PetscFunctionReturn(PETSC_SUCCESS);
3557: }
3559: /*@C
3560: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3562: Not Collective
3564: Input Parameters:
3565: + section - the section
3566: . field - the field number
3567: . numPoints - the number of points
3568: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3569: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3570: context, see `DMPlexGetConeOrientation()`).
3571: . perms - The permutations for the given orientations: set to NULL at conclusion
3572: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3574: Level: developer
3576: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3577: @*/
3578: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3579: {
3580: PetscFunctionBegin;
3582: 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);
3583: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3584: PetscFunctionReturn(PETSC_SUCCESS);
3585: }
3587: /*@
3588: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3590: Not Collective
3592: Input Parameter:
3593: . sym - the `PetscSectionSym`
3595: Output Parameter:
3596: . nsym - the equivalent symmetries
3598: Level: developer
3600: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3601: @*/
3602: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3603: {
3604: PetscFunctionBegin;
3607: PetscTryTypeMethod(sym, copy, nsym);
3608: PetscFunctionReturn(PETSC_SUCCESS);
3609: }
3611: /*@
3612: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3614: Collective
3616: Input Parameters:
3617: + sym - the `PetscSectionSym`
3618: - migrationSF - the distribution map from roots to leaves
3620: Output Parameter:
3621: . dsym - the redistributed symmetries
3623: Level: developer
3625: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3626: @*/
3627: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3628: {
3629: PetscFunctionBegin;
3633: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3634: PetscFunctionReturn(PETSC_SUCCESS);
3635: }
3637: /*@
3638: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3640: Not Collective
3642: Input Parameter:
3643: . s - the global `PetscSection`
3645: Output Parameter:
3646: . flg - the flag
3648: Level: developer
3650: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3651: @*/
3652: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3653: {
3654: PetscFunctionBegin;
3656: *flg = s->useFieldOff;
3657: PetscFunctionReturn(PETSC_SUCCESS);
3658: }
3660: /*@
3661: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3663: Not Collective
3665: Input Parameters:
3666: + s - the global `PetscSection`
3667: - flg - the flag
3669: Level: developer
3671: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3672: @*/
3673: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3674: {
3675: PetscFunctionBegin;
3677: s->useFieldOff = flg;
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: #define PetscSectionExpandPoints_Loop(TYPE) \
3682: { \
3683: PetscInt i, n, o0, o1, size; \
3684: TYPE *a0 = (TYPE *)origArray, *a1; \
3685: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3686: PetscCall(PetscMalloc1(size, &a1)); \
3687: for (i = 0; i < npoints; i++) { \
3688: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3689: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3690: PetscCall(PetscSectionGetDof(s, i, &n)); \
3691: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3692: } \
3693: *newArray = (void *)a1; \
3694: }
3696: /*@
3697: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3699: Not Collective
3701: Input Parameters:
3702: + origSection - the `PetscSection` describing the layout of the array
3703: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3704: . origArray - the array; its size must be equal to the storage size of `origSection`
3705: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
3707: Output Parameters:
3708: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3709: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
3711: Level: developer
3713: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3714: @*/
3715: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3716: {
3717: PetscSection s;
3718: const PetscInt *points_;
3719: PetscInt i, n, npoints, pStart, pEnd;
3720: PetscMPIInt unitsize;
3722: PetscFunctionBegin;
3728: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3729: PetscCall(ISGetLocalSize(points, &npoints));
3730: PetscCall(ISGetIndices(points, &points_));
3731: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3732: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3733: PetscCall(PetscSectionSetChart(s, 0, npoints));
3734: for (i = 0; i < npoints; i++) {
3735: 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);
3736: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3737: PetscCall(PetscSectionSetDof(s, i, n));
3738: }
3739: PetscCall(PetscSectionSetUp(s));
3740: if (newArray) {
3741: if (dataType == MPIU_INT) {
3742: PetscSectionExpandPoints_Loop(PetscInt);
3743: } else if (dataType == MPIU_SCALAR) {
3744: PetscSectionExpandPoints_Loop(PetscScalar);
3745: } else if (dataType == MPIU_REAL) {
3746: PetscSectionExpandPoints_Loop(PetscReal);
3747: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3748: }
3749: if (newSection) {
3750: *newSection = s;
3751: } else {
3752: PetscCall(PetscSectionDestroy(&s));
3753: }
3754: PetscCall(ISRestoreIndices(points, &points_));
3755: PetscFunctionReturn(PETSC_SUCCESS);
3756: }