Actual source code: section.c

petsc-main 2021-04-20
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

  5: #include <petsc/private/sectionimpl.h>
  6: #include <petscsection.h>
  7: #include <petscsf.h>
  8: #include <petscviewer.h>

 10: PetscClassId PETSC_SECTION_CLASSID;

 12: /*@
 13:   PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.

 15:   Collective

 17:   Input Parameters:
 18: + comm - the MPI communicator
 19: - s    - pointer to the section

 21:   Level: beginner

 23:   Notes:
 24:   Typical calling sequence
 25: $       PetscSectionCreate(MPI_Comm,PetscSection *);
 26: $       PetscSectionSetNumFields(PetscSection, numFields);
 27: $       PetscSectionSetChart(PetscSection,low,high);
 28: $       PetscSectionSetDof(PetscSection,point,numdof);
 29: $       PetscSectionSetUp(PetscSection);
 30: $       PetscSectionGetOffset(PetscSection,point,PetscInt *);
 31: $       PetscSectionDestroy(PetscSection);

 33:   The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
 34:   recommended they not be used in user codes unless you really gain something in their use.

 36: .seealso: PetscSection, PetscSectionDestroy()
 37: @*/
 38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
 39: {

 44:   ISInitializePackage();

 46:   PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);

 48:   (*s)->pStart              = -1;
 49:   (*s)->pEnd                = -1;
 50:   (*s)->perm                = NULL;
 51:   (*s)->pointMajor          = PETSC_TRUE;
 52:   (*s)->includesConstraints = PETSC_TRUE;
 53:   (*s)->maxDof              = 0;
 54:   (*s)->atlasDof            = NULL;
 55:   (*s)->atlasOff            = NULL;
 56:   (*s)->bc                  = NULL;
 57:   (*s)->bcIndices           = NULL;
 58:   (*s)->setup               = PETSC_FALSE;
 59:   (*s)->numFields           = 0;
 60:   (*s)->fieldNames          = NULL;
 61:   (*s)->field               = NULL;
 62:   (*s)->useFieldOff         = PETSC_FALSE;
 63:   (*s)->compNames           = NULL;
 64:   (*s)->clObj               = NULL;
 65:   (*s)->clHash              = NULL;
 66:   (*s)->clSection           = NULL;
 67:   (*s)->clPoints            = NULL;
 68:   return(0);
 69: }

 71: /*@
 72:   PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection

 74:   Collective

 76:   Input Parameter:
 77: . section - the PetscSection

 79:   Output Parameter:
 80: . newSection - the copy

 82:   Level: intermediate

 84: .seealso: 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;
 91:   PetscErrorCode  ierr;

 96:   PetscSectionReset(newSection);
 97:   PetscSectionGetNumFields(section, &numFields);
 98:   if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
 99:   for (f = 0; f < numFields; ++f) {
100:     const char *fieldName = NULL, *compName = NULL;
101:     PetscInt   numComp = 0;

103:     PetscSectionGetFieldName(section, f, &fieldName);
104:     PetscSectionSetFieldName(newSection, f, fieldName);
105:     PetscSectionGetFieldComponents(section, f, &numComp);
106:     PetscSectionSetFieldComponents(newSection, f, numComp);
107:     for (c = 0; c < numComp; ++c) {
108:       PetscSectionGetComponentName(section, f, c, &compName);
109:       PetscSectionSetComponentName(newSection, f, c, compName);
110:     }
111:     PetscSectionGetFieldSym(section, f, &sym);
112:     PetscSectionSetFieldSym(newSection, f, sym);
113:   }
114:   PetscSectionGetChart(section, &pStart, &pEnd);
115:   PetscSectionSetChart(newSection, pStart, pEnd);
116:   PetscSectionGetPermutation(section, &perm);
117:   PetscSectionSetPermutation(newSection, perm);
118:   PetscSectionGetSym(section, &sym);
119:   PetscSectionSetSym(newSection, sym);
120:   for (p = pStart; p < pEnd; ++p) {
121:     PetscInt dof, cdof, fcdof = 0;

123:     PetscSectionGetDof(section, p, &dof);
124:     PetscSectionSetDof(newSection, p, dof);
125:     PetscSectionGetConstraintDof(section, p, &cdof);
126:     if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
127:     for (f = 0; f < numFields; ++f) {
128:       PetscSectionGetFieldDof(section, p, f, &dof);
129:       PetscSectionSetFieldDof(newSection, p, f, dof);
130:       if (cdof) {
131:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
132:         if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
133:       }
134:     }
135:   }
136:   PetscSectionSetUp(newSection);
137:   for (p = pStart; p < pEnd; ++p) {
138:     PetscInt       off, cdof, fcdof = 0;
139:     const PetscInt *cInd;

141:     /* Must set offsets in case they do not agree with the prefix sums */
142:     PetscSectionGetOffset(section, p, &off);
143:     PetscSectionSetOffset(newSection, p, off);
144:     PetscSectionGetConstraintDof(section, p, &cdof);
145:     if (cdof) {
146:       PetscSectionGetConstraintIndices(section, p, &cInd);
147:       PetscSectionSetConstraintIndices(newSection, p, cInd);
148:       for (f = 0; f < numFields; ++f) {
149:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
150:         if (fcdof) {
151:           PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
152:           PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
153:         }
154:       }
155:     }
156:   }
157:   return(0);
158: }

160: /*@
161:   PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection

163:   Collective

165:   Input Parameter:
166: . section - the PetscSection

168:   Output Parameter:
169: . newSection - the copy

171:   Level: beginner

173: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
174: @*/
175: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
176: {

182:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
183:   PetscSectionCopy(section, *newSection);
184:   return(0);
185: }

187: /*@
188:   PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database

190:   Collective on PetscSection

192:   Input Parameter:
193: . section - the PetscSection

195:   Options Database:
196: . -petscsection_point_major the dof order

198:   Level: intermediate

200: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
201: @*/
202: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
203: {

208:   PetscObjectOptionsBegin((PetscObject) s);
209:   PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
210:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
211:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
212:   PetscOptionsEnd();
213:   PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
214:   return(0);
215: }

217: /*@
218:   PetscSectionCompare - Compares two sections

220:   Collective on PetscSection

222:   Input Parameters:
223: + s1 - the first PetscSection
224: - s2 - the second PetscSection

226:   Output Parameter:
227: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise

229:   Level: intermediate

231:   Notes:
232:   Field names are disregarded.

234: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
235: @*/
236: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
237: {
238:   PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
239:   const PetscInt *idx1, *idx2;
240:   IS perm1, perm2;
241:   PetscBool flg;
242:   PetscMPIInt mflg;

249:   flg = PETSC_FALSE;

251:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
252:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
253:     *congruent = PETSC_FALSE;
254:     return(0);
255:   }

257:   PetscSectionGetChart(s1, &pStart, &pEnd);
258:   PetscSectionGetChart(s2, &n1, &n2);
259:   if (pStart != n1 || pEnd != n2) goto not_congruent;

261:   PetscSectionGetPermutation(s1, &perm1);
262:   PetscSectionGetPermutation(s2, &perm2);
263:   if (perm1 && perm2) {
264:     ISEqual(perm1, perm2, congruent);
265:     if (!(*congruent)) goto not_congruent;
266:   } else if (perm1 != perm2) goto not_congruent;

268:   for (p = pStart; p < pEnd; ++p) {
269:     PetscSectionGetOffset(s1, p, &n1);
270:     PetscSectionGetOffset(s2, p, &n2);
271:     if (n1 != n2) goto not_congruent;

273:     PetscSectionGetDof(s1, p, &n1);
274:     PetscSectionGetDof(s2, p, &n2);
275:     if (n1 != n2) goto not_congruent;

277:     PetscSectionGetConstraintDof(s1, p, &ncdof);
278:     PetscSectionGetConstraintDof(s2, p, &n2);
279:     if (ncdof != n2) goto not_congruent;

281:     PetscSectionGetConstraintIndices(s1, p, &idx1);
282:     PetscSectionGetConstraintIndices(s2, p, &idx2);
283:     PetscArraycmp(idx1, idx2, ncdof, congruent);
284:     if (!(*congruent)) goto not_congruent;
285:   }

287:   PetscSectionGetNumFields(s1, &nfields);
288:   PetscSectionGetNumFields(s2, &n2);
289:   if (nfields != n2) goto not_congruent;

291:   for (f = 0; f < nfields; ++f) {
292:     PetscSectionGetFieldComponents(s1, f, &n1);
293:     PetscSectionGetFieldComponents(s2, f, &n2);
294:     if (n1 != n2) goto not_congruent;

296:     for (p = pStart; p < pEnd; ++p) {
297:       PetscSectionGetFieldOffset(s1, p, f, &n1);
298:       PetscSectionGetFieldOffset(s2, p, f, &n2);
299:       if (n1 != n2) goto not_congruent;

301:       PetscSectionGetFieldDof(s1, p, f, &n1);
302:       PetscSectionGetFieldDof(s2, p, f, &n2);
303:       if (n1 != n2) goto not_congruent;

305:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
306:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
307:       if (nfcdof != n2) goto not_congruent;

309:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
310:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
311:       PetscArraycmp(idx1, idx2, nfcdof, congruent);
312:       if (!(*congruent)) goto not_congruent;
313:     }
314:   }

316:   flg = PETSC_TRUE;
317: not_congruent:
318:   MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
319:   return(0);
320: }

322: /*@
323:   PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.

325:   Not collective

327:   Input Parameter:
328: . s - the PetscSection

330:   Output Parameter:
331: . numFields - the number of fields defined, or 0 if none were defined

333:   Level: intermediate

335: .seealso: PetscSectionSetNumFields()
336: @*/
337: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
338: {
342:   *numFields = s->numFields;
343:   return(0);
344: }

346: /*@
347:   PetscSectionSetNumFields - Sets the number of fields.

349:   Not collective

351:   Input Parameters:
352: + s - the PetscSection
353: - numFields - the number of fields

355:   Level: intermediate

357: .seealso: PetscSectionGetNumFields()
358: @*/
359: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
360: {
361:   PetscInt       f;

366:   if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
367:   PetscSectionReset(s);

369:   s->numFields = numFields;
370:   PetscMalloc1(s->numFields, &s->numFieldComponents);
371:   PetscMalloc1(s->numFields, &s->fieldNames);
372:   PetscMalloc1(s->numFields, &s->compNames);
373:   PetscMalloc1(s->numFields, &s->field);
374:   for (f = 0; f < s->numFields; ++f) {
375:     char name[64];

377:     s->numFieldComponents[f] = 1;

379:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
380:     PetscSNPrintf(name, 64, "Field_%D", f);
381:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
382:     PetscSNPrintf(name, 64, "Component_0");
383:     PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
384:     PetscStrallocpy(name, (char **) &s->compNames[f][0]);
385:   }
386:   return(0);
387: }

389: /*@C
390:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

392:   Not Collective

394:   Input Parameters:
395: + s     - the PetscSection
396: - field - the field number

398:   Output Parameter:
399: . fieldName - the field name

401:   Level: intermediate

403: .seealso: PetscSectionSetFieldName()
404: @*/
405: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
406: {
410:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
411:   *fieldName = s->fieldNames[field];
412:   return(0);
413: }

415: /*@C
416:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

418:   Not Collective

420:   Input Parameters:
421: + s     - the PetscSection
422: . field - the field number
423: - fieldName - the field name

425:   Level: intermediate

427: .seealso: PetscSectionGetFieldName()
428: @*/
429: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
430: {

436:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
437:   PetscFree(s->fieldNames[field]);
438:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
439:   return(0);
440: }

442: /*@C
443:   PetscSectionGetComponentName - Gets the name of a field component in the PetscSection

445:   Not Collective

447:   Input Parameters:
448: + s     - the PetscSection
449: . field - the field number
450: . comp  - the component number
451: - compName - the component name

453:   Level: intermediate

455: .seealso: PetscSectionSetComponentName()
456: @*/
457: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
458: {
462:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
463:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
464:   *compName = s->compNames[field][comp];
465:   return(0);
466: }

468: /*@C
469:   PetscSectionSetComponentName - Sets the name of a field component in the PetscSection

471:   Not Collective

473:   Input Parameters:
474: + s     - the PetscSection
475: . field - the field number
476: . comp  - the component number
477: - compName - the component name

479:   Level: intermediate

481: .seealso: PetscSectionGetComponentName()
482: @*/
483: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
484: {

490:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
491:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
492:   PetscFree(s->compNames[field][comp]);
493:   PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
494:   return(0);
495: }

497: /*@
498:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

500:   Not collective

502:   Input Parameters:
503: + s - the PetscSection
504: - field - the field number

506:   Output Parameter:
507: . numComp - the number of field components

509:   Level: intermediate

511: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
512: @*/
513: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
514: {
518:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
519:   *numComp = s->numFieldComponents[field];
520:   return(0);
521: }

523: /*@
524:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

526:   Not collective

528:   Input Parameters:
529: + s - the PetscSection
530: . field - the field number
531: - numComp - the number of field components

533:   Level: intermediate

535: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
536: @*/
537: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
538: {
540:   PetscInt c;
541:   char name[64];

545:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
546:   if (s->compNames) {
547:     for (c = 0; c < s->numFieldComponents[field]; ++c) {
548:       PetscFree(s->compNames[field][c]);
549:     }
550:     PetscFree(s->compNames[field]);
551:   }

553:   s->numFieldComponents[field] = numComp;
554:   if (numComp) {
555:     PetscMalloc1(numComp, (char ***) &s->compNames[field]);
556:     for (c = 0; c < numComp; ++c) {
557:       PetscSNPrintf(name, 64, "%D", c);
558:       PetscStrallocpy(name, (char **) &s->compNames[field][c]);
559:     }
560:   }
561:   return(0);
562: }

564: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
565: {

569:   if (!s->bc) {
570:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
571:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
572:   }
573:   return(0);
574: }

576: /*@
577:   PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.

579:   Not collective

581:   Input Parameter:
582: . s - the PetscSection

584:   Output Parameters:
585: + pStart - the first point
586: - pEnd - one past the last point

588:   Level: intermediate

590: .seealso: PetscSectionSetChart(), PetscSectionCreate()
591: @*/
592: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
593: {
596:   if (pStart) *pStart = s->pStart;
597:   if (pEnd)   *pEnd   = s->pEnd;
598:   return(0);
599: }

601: /*@
602:   PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.

604:   Not collective

606:   Input Parameters:
607: + s - the PetscSection
608: . pStart - the first point
609: - pEnd - one past the last point

611:   Level: intermediate

613: .seealso: PetscSectionGetChart(), PetscSectionCreate()
614: @*/
615: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
616: {
617:   PetscInt       f;

622:   /* Cannot Reset() because it destroys field information */
623:   s->setup = PETSC_FALSE;
624:   PetscSectionDestroy(&s->bc);
625:   PetscFree(s->bcIndices);
626:   PetscFree2(s->atlasDof, s->atlasOff);

628:   s->pStart = pStart;
629:   s->pEnd   = pEnd;
630:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
631:   PetscArrayzero(s->atlasDof, pEnd - pStart);
632:   for (f = 0; f < s->numFields; ++f) {
633:     PetscSectionSetChart(s->field[f], pStart, pEnd);
634:   }
635:   return(0);
636: }

638: /*@
639:   PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL

641:   Not collective

643:   Input Parameter:
644: . s - the PetscSection

646:   Output Parameters:
647: . perm - The permutation as an IS

649:   Level: intermediate

651: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
652: @*/
653: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
654: {
658:   return(0);
659: }

661: /*@
662:   PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)

664:   Not collective

666:   Input Parameters:
667: + s - the PetscSection
668: - perm - the permutation of points

670:   Level: intermediate

672: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
673: @*/
674: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
675: {

681:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
682:   if (s->perm != perm) {
683:     ISDestroy(&s->perm);
684:     if (perm) {
685:       s->perm = perm;
686:       PetscObjectReference((PetscObject) s->perm);
687:     }
688:   }
689:   return(0);
690: }

692: /*@
693:   PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major

695:   Not collective

697:   Input Parameter:
698: . s - the PetscSection

700:   Output Parameter:
701: . pm - the flag for point major ordering

703:   Level: intermediate

705: .seealso: PetscSectionSetPointMajor()
706: @*/
707: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
708: {
712:   *pm = s->pointMajor;
713:   return(0);
714: }

716: /*@
717:   PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major

719:   Not collective

721:   Input Parameters:
722: + s  - the PetscSection
723: - pm - the flag for point major ordering

725:   Not collective

727:   Level: intermediate

729: .seealso: PetscSectionGetPointMajor()
730: @*/
731: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
732: {
735:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
736:   s->pointMajor = pm;
737:   return(0);
738: }

740: /*@
741:   PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets

743:   Not collective

745:   Input Parameter:
746: . s - the PetscSection

748:   Output Parameter:
749: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets

751:   Level: intermediate

753: .seealso: PetscSectionSetIncludesConstraints()
754: @*/
755: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
756: {
760:   *includesConstraints = s->includesConstraints;
761:   return(0);
762: }

764: /*@
765:   PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets

767:   Not collective

769:   Input Parameters:
770: + s  - the PetscSection
771: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets

773:   Not collective

775:   Level: intermediate

777: .seealso: PetscSectionGetIncludesConstraints()
778: @*/
779: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
780: {
783:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
784:   s->includesConstraints = includesConstraints;
785:   return(0);
786: }

788: /*@
789:   PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.

791:   Not collective

793:   Input Parameters:
794: + s - the PetscSection
795: - point - the point

797:   Output Parameter:
798: . numDof - the number of dof

800:   Level: intermediate

802: .seealso: PetscSectionSetDof(), PetscSectionCreate()
803: @*/
804: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
805: {
809:   if (PetscDefined(USE_DEBUG)) {
810:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
811:   }
812:   *numDof = s->atlasDof[point - s->pStart];
813:   return(0);
814: }

816: /*@
817:   PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.

819:   Not collective

821:   Input Parameters:
822: + s - the PetscSection
823: . point - the point
824: - numDof - the number of dof

826:   Level: intermediate

828: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
829: @*/
830: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
831: {
834:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
835:   s->atlasDof[point - s->pStart] = numDof;
836:   return(0);
837: }

839: /*@
840:   PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.

842:   Not collective

844:   Input Parameters:
845: + s - the PetscSection
846: . point - the point
847: - numDof - the number of additional dof

849:   Level: intermediate

851: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
852: @*/
853: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
854: {
857:   if (PetscDefined(USE_DEBUG)) {
858:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
859:   }
860:   s->atlasDof[point - s->pStart] += numDof;
861:   return(0);
862: }

864: /*@
865:   PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.

867:   Not collective

869:   Input Parameters:
870: + s - the PetscSection
871: . point - the point
872: - field - the field

874:   Output Parameter:
875: . numDof - the number of dof

877:   Level: intermediate

879: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
880: @*/
881: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
882: {

887:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
888:   PetscSectionGetDof(s->field[field], point, numDof);
889:   return(0);
890: }

892: /*@
893:   PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.

895:   Not collective

897:   Input Parameters:
898: + s - the PetscSection
899: . point - the point
900: . field - the field
901: - numDof - the number of dof

903:   Level: intermediate

905: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
906: @*/
907: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
908: {

913:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
914:   PetscSectionSetDof(s->field[field], point, numDof);
915:   return(0);
916: }

918: /*@
919:   PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.

921:   Not collective

923:   Input Parameters:
924: + s - the PetscSection
925: . point - the point
926: . field - the field
927: - numDof - the number of dof

929:   Level: intermediate

931: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
932: @*/
933: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
934: {

939:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
940:   PetscSectionAddDof(s->field[field], point, numDof);
941:   return(0);
942: }

944: /*@
945:   PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.

947:   Not collective

949:   Input Parameters:
950: + s - the PetscSection
951: - point - the point

953:   Output Parameter:
954: . numDof - the number of dof which are fixed by constraints

956:   Level: intermediate

958: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
959: @*/
960: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
961: {

967:   if (s->bc) {
968:     PetscSectionGetDof(s->bc, point, numDof);
969:   } else *numDof = 0;
970:   return(0);
971: }

973: /*@
974:   PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.

976:   Not collective

978:   Input Parameters:
979: + s - the PetscSection
980: . point - the point
981: - numDof - the number of dof which are fixed by constraints

983:   Level: intermediate

985: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
986: @*/
987: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
988: {

993:   if (numDof) {
994:     PetscSectionCheckConstraints_Static(s);
995:     PetscSectionSetDof(s->bc, point, numDof);
996:   }
997:   return(0);
998: }

1000: /*@
1001:   PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.

1003:   Not collective

1005:   Input Parameters:
1006: + s - the PetscSection
1007: . point - the point
1008: - numDof - the number of additional dof which are fixed by constraints

1010:   Level: intermediate

1012: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
1013: @*/
1014: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1015: {

1020:   if (numDof) {
1021:     PetscSectionCheckConstraints_Static(s);
1022:     PetscSectionAddDof(s->bc, point, numDof);
1023:   }
1024:   return(0);
1025: }

1027: /*@
1028:   PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.

1030:   Not collective

1032:   Input Parameters:
1033: + s - the PetscSection
1034: . point - the point
1035: - field - the field

1037:   Output Parameter:
1038: . numDof - the number of dof which are fixed by constraints

1040:   Level: intermediate

1042: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
1043: @*/
1044: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1045: {

1051:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1052:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
1053:   return(0);
1054: }

1056: /*@
1057:   PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.

1059:   Not collective

1061:   Input Parameters:
1062: + s - the PetscSection
1063: . point - the point
1064: . field - the field
1065: - numDof - the number of dof which are fixed by constraints

1067:   Level: intermediate

1069: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1070: @*/
1071: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1072: {

1077:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1078:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
1079:   return(0);
1080: }

1082: /*@
1083:   PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.

1085:   Not collective

1087:   Input Parameters:
1088: + s - the PetscSection
1089: . point - the point
1090: . field - the field
1091: - numDof - the number of additional dof which are fixed by constraints

1093:   Level: intermediate

1095: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1096: @*/
1097: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1098: {

1103:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1104:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
1105:   return(0);
1106: }

1108: /*@
1109:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1111:   Not collective

1113:   Input Parameter:
1114: . s - the PetscSection

1116:   Level: advanced

1118: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1119: @*/
1120: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1121: {

1126:   if (s->bc) {
1127:     const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;

1129:     PetscSectionSetUp(s->bc);
1130:     PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1131:   }
1132:   return(0);
1133: }

1135: /*@
1136:   PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.

1138:   Not collective

1140:   Input Parameter:
1141: . s - the PetscSection

1143:   Level: intermediate

1145: .seealso: PetscSectionCreate()
1146: @*/
1147: PetscErrorCode PetscSectionSetUp(PetscSection s)
1148: {
1149:   const PetscInt *pind   = NULL;
1150:   PetscInt        offset = 0, foff, p, f;
1151:   PetscErrorCode  ierr;

1155:   if (s->setup) return(0);
1156:   s->setup = PETSC_TRUE;
1157:   /* Set offsets and field offsets for all points */
1158:   /*   Assume that all fields have the same chart */
1159:   if (!s->includesConstraints) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1160:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1161:   if (s->pointMajor) {
1162:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1163:       const PetscInt q = pind ? pind[p] : p;

1165:       /* Set point offset */
1166:       s->atlasOff[q] = offset;
1167:       offset        += s->atlasDof[q];
1168:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
1169:       /* Set field offset */
1170:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1171:         PetscSection sf = s->field[f];

1173:         sf->atlasOff[q] = foff;
1174:         foff += sf->atlasDof[q];
1175:       }
1176:     }
1177:   } else {
1178:     /* Set field offsets for all points */
1179:     for (f = 0; f < s->numFields; ++f) {
1180:       PetscSection sf = s->field[f];

1182:       for (p = 0; p < s->pEnd - s->pStart; ++p) {
1183:         const PetscInt q = pind ? pind[p] : p;

1185:         sf->atlasOff[q] = offset;
1186:         offset += sf->atlasDof[q];
1187:       }
1188:     }
1189:     /* Disable point offsets since these are unused */
1190:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1191:       s->atlasOff[p] = -1;
1192:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
1193:     }
1194:   }
1195:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1196:   /* Setup BC sections */
1197:   PetscSectionSetUpBC(s);
1198:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1199:   return(0);
1200: }

1202: /*@
1203:   PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart

1205:   Not collective

1207:   Input Parameters:
1208: . s - the PetscSection

1210:   Output Parameter:
1211: . maxDof - the maximum dof

1213:   Level: intermediate

1215: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1216: @*/
1217: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1218: {
1222:   *maxDof = s->maxDof;
1223:   return(0);
1224: }

1226: /*@
1227:   PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.

1229:   Not collective

1231:   Input Parameter:
1232: . s - the PetscSection

1234:   Output Parameter:
1235: . size - the size of an array which can hold all the dofs

1237:   Level: intermediate

1239: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1240: @*/
1241: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1242: {
1243:   PetscInt p, n = 0;

1248:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1249:   *size = n;
1250:   return(0);
1251: }

1253: /*@
1254:   PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.

1256:   Not collective

1258:   Input Parameter:
1259: . s - the PetscSection

1261:   Output Parameter:
1262: . size - the size of an array which can hold all unconstrained dofs

1264:   Level: intermediate

1266: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1267: @*/
1268: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1269: {
1270:   PetscInt p, n = 0;

1275:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1276:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1277:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1278:   }
1279:   *size = n;
1280:   return(0);
1281: }

1283: /*@
1284:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1285:   the local section and an SF describing the section point overlap.

1287:   Input Parameters:
1288: + s - The PetscSection for the local field layout
1289: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1290: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1291: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points

1293:   Output Parameter:
1294: . gsection - The PetscSection for the global field layout

1296:   Note: This gives negative sizes and offsets to points not owned by this process

1298:   Level: intermediate

1300: .seealso: PetscSectionCreate()
1301: @*/
1302: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1303: {
1304:   PetscSection    gs;
1305:   const PetscInt *pind = NULL;
1306:   PetscInt       *recv = NULL, *neg = NULL;
1307:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1308:   PetscErrorCode  ierr;

1316:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1317:   PetscSectionGetChart(s, &pStart, &pEnd);
1318:   PetscSectionSetChart(gs, pStart, pEnd);
1319:   gs->includesConstraints = includeConstraints;
1320:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1321:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1322:   /* Must allocate for all points visible to SF, which may be more than this section */
1323:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1324:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1325:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1326:     if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1327:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1328:     PetscArrayzero(neg,nroots);
1329:   }
1330:   /* Mark all local points with negative dof */
1331:   for (p = pStart; p < pEnd; ++p) {
1332:     PetscSectionGetDof(s, p, &dof);
1333:     PetscSectionSetDof(gs, p, dof);
1334:     PetscSectionGetConstraintDof(s, p, &cdof);
1335:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1336:     if (neg) neg[p] = -(dof+1);
1337:   }
1338:   PetscSectionSetUpBC(gs);
1339:   if (gs->bcIndices) {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]);}
1340:   if (nroots >= 0) {
1341:     PetscArrayzero(recv,nlocal);
1342:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1343:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1344:     for (p = pStart; p < pEnd; ++p) {
1345:       if (recv[p] < 0) {
1346:         gs->atlasDof[p-pStart] = recv[p];
1347:         PetscSectionGetDof(s, p, &dof);
1348:         if (-(recv[p]+1) != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %D for point %D is not the unconstrained %D", -(recv[p]+1), p, dof);
1349:       }
1350:     }
1351:   }
1352:   /* Calculate new sizes, get process offset, and calculate point offsets */
1353:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1354:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1355:     const PetscInt q = pind ? pind[p] : p;

1357:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1358:     gs->atlasOff[q] = off;
1359:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1360:   }
1361:   if (!localOffsets) {
1362:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1363:     globalOff -= off;
1364:   }
1365:   for (p = pStart, off = 0; p < pEnd; ++p) {
1366:     gs->atlasOff[p-pStart] += globalOff;
1367:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1368:   }
1369:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1370:   /* Put in negative offsets for ghost points */
1371:   if (nroots >= 0) {
1372:     PetscArrayzero(recv,nlocal);
1373:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1374:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1375:     for (p = pStart; p < pEnd; ++p) {
1376:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1377:     }
1378:   }
1379:   PetscFree2(neg,recv);
1380:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1381:   *gsection = gs;
1382:   return(0);
1383: }

1385: /*@
1386:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1387:   the local section and an SF describing the section point overlap.

1389:   Input Parameters:
1390:   + s - The PetscSection for the local field layout
1391:   . sf - The SF describing parallel layout of the section points
1392:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1393:   . numExcludes - The number of exclusion ranges
1394:   - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs

1396:   Output Parameter:
1397:   . gsection - The PetscSection for the global field layout

1399:   Note: This gives negative sizes and offsets to points not owned by this process

1401:   Level: advanced

1403: .seealso: PetscSectionCreate()
1404: @*/
1405: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1406: {
1407:   const PetscInt *pind = NULL;
1408:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1409:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1410:   PetscErrorCode  ierr;

1416:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1417:   PetscSectionGetChart(s, &pStart, &pEnd);
1418:   PetscSectionSetChart(*gsection, pStart, pEnd);
1419:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1420:   if (nroots >= 0) {
1421:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1422:     PetscCalloc1(nroots, &neg);
1423:     if (nroots > pEnd-pStart) {
1424:       PetscCalloc1(nroots, &tmpOff);
1425:     } else {
1426:       tmpOff = &(*gsection)->atlasDof[-pStart];
1427:     }
1428:   }
1429:   /* Mark ghost points with negative dof */
1430:   for (p = pStart; p < pEnd; ++p) {
1431:     for (e = 0; e < numExcludes; ++e) {
1432:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1433:         PetscSectionSetDof(*gsection, p, 0);
1434:         break;
1435:       }
1436:     }
1437:     if (e < numExcludes) continue;
1438:     PetscSectionGetDof(s, p, &dof);
1439:     PetscSectionSetDof(*gsection, p, dof);
1440:     PetscSectionGetConstraintDof(s, p, &cdof);
1441:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1442:     if (neg) neg[p] = -(dof+1);
1443:   }
1444:   PetscSectionSetUpBC(*gsection);
1445:   if (nroots >= 0) {
1446:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1447:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1448:     if (nroots > pEnd - pStart) {
1449:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1450:     }
1451:   }
1452:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1453:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1454:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1455:     const PetscInt q = pind ? pind[p] : p;

1457:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1458:     (*gsection)->atlasOff[q] = off;
1459:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1460:   }
1461:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1462:   globalOff -= off;
1463:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1464:     (*gsection)->atlasOff[p] += globalOff;
1465:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1466:   }
1467:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1468:   /* Put in negative offsets for ghost points */
1469:   if (nroots >= 0) {
1470:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1471:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1472:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1473:     if (nroots > pEnd - pStart) {
1474:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1475:     }
1476:   }
1477:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1478:   PetscFree(neg);
1479:   return(0);
1480: }

1482: /*@
1483:   PetscSectionGetPointLayout - Get the PetscLayout associated with the section points

1485:   Collective on comm

1487:   Input Parameters:
1488: + comm - The MPI_Comm
1489: - s    - The PetscSection

1491:   Output Parameter:
1492: . layout - The point layout for the section

1494:   Note: This is usually caleld for the default global section.

1496:   Level: advanced

1498: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1499: @*/
1500: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1501: {
1502:   PetscInt       pStart, pEnd, p, localSize = 0;

1506:   PetscSectionGetChart(s, &pStart, &pEnd);
1507:   for (p = pStart; p < pEnd; ++p) {
1508:     PetscInt dof;

1510:     PetscSectionGetDof(s, p, &dof);
1511:     if (dof > 0) ++localSize;
1512:   }
1513:   PetscLayoutCreate(comm, layout);
1514:   PetscLayoutSetLocalSize(*layout, localSize);
1515:   PetscLayoutSetBlockSize(*layout, 1);
1516:   PetscLayoutSetUp(*layout);
1517:   return(0);
1518: }

1520: /*@
1521:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

1523:   Collective on comm

1525:   Input Parameters:
1526: + comm - The MPI_Comm
1527: - s    - The PetscSection

1529:   Output Parameter:
1530: . layout - The dof layout for the section

1532:   Note: This is usually called for the default global section.

1534:   Level: advanced

1536: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1537: @*/
1538: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1539: {
1540:   PetscInt       pStart, pEnd, p, localSize = 0;

1546:   PetscSectionGetChart(s, &pStart, &pEnd);
1547:   for (p = pStart; p < pEnd; ++p) {
1548:     PetscInt dof,cdof;

1550:     PetscSectionGetDof(s, p, &dof);
1551:     PetscSectionGetConstraintDof(s, p, &cdof);
1552:     if (dof-cdof > 0) localSize += dof-cdof;
1553:   }
1554:   PetscLayoutCreate(comm, layout);
1555:   PetscLayoutSetLocalSize(*layout, localSize);
1556:   PetscLayoutSetBlockSize(*layout, 1);
1557:   PetscLayoutSetUp(*layout);
1558:   return(0);
1559: }

1561: /*@
1562:   PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.

1564:   Not collective

1566:   Input Parameters:
1567: + s - the PetscSection
1568: - point - the point

1570:   Output Parameter:
1571: . offset - the offset

1573:   Level: intermediate

1575: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1576: @*/
1577: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1578: {
1582:   if (PetscDefined(USE_DEBUG)) {
1583:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
1584:   }
1585:   *offset = s->atlasOff[point - s->pStart];
1586:   return(0);
1587: }

1589: /*@
1590:   PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.

1592:   Not collective

1594:   Input Parameters:
1595: + s - the PetscSection
1596: . point - the point
1597: - offset - the offset

1599:   Note: The user usually does not call this function, but uses PetscSectionSetUp()

1601:   Level: intermediate

1603: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1604: @*/
1605: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1606: {
1609:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
1610:   s->atlasOff[point - s->pStart] = offset;
1611:   return(0);
1612: }

1614: /*@
1615:   PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.

1617:   Not collective

1619:   Input Parameters:
1620: + s - the PetscSection
1621: . point - the point
1622: - field - the field

1624:   Output Parameter:
1625: . offset - the offset

1627:   Level: intermediate

1629: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1630: @*/
1631: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1632: {

1638:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1639:   PetscSectionGetOffset(s->field[field], point, offset);
1640:   return(0);
1641: }

1643: /*@
1644:   PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.

1646:   Not collective

1648:   Input Parameters:
1649: + s - the PetscSection
1650: . point - the point
1651: . field - the field
1652: - offset - the offset

1654:   Note: The user usually does not call this function, but uses PetscSectionSetUp()

1656:   Level: intermediate

1658: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1659: @*/
1660: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1661: {

1666:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1667:   PetscSectionSetOffset(s->field[field], point, offset);
1668:   return(0);
1669: }

1671: /*@
1672:   PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.

1674:   Not collective

1676:   Input Parameters:
1677: + s - the PetscSection
1678: . point - the point
1679: - field - the field

1681:   Output Parameter:
1682: . offset - the offset

1684:   Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1685:         this point, what number is the first dof with this field.

1687:   Level: advanced

1689: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1690: @*/
1691: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1692: {
1693:   PetscInt       off, foff;

1699:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
1700:   PetscSectionGetOffset(s, point, &off);
1701:   PetscSectionGetOffset(s->field[field], point, &foff);
1702:   *offset = foff - off;
1703:   return(0);
1704: }

1706: /*@
1707:   PetscSectionGetOffsetRange - Return the full range of offsets [start, end)

1709:   Not collective

1711:   Input Parameter:
1712: . s - the PetscSection

1714:   Output Parameters:
1715: + start - the minimum offset
1716: - end   - one more than the maximum offset

1718:   Level: intermediate

1720: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1721: @*/
1722: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1723: {
1724:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1729:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1730:   PetscSectionGetChart(s, &pStart, &pEnd);
1731:   for (p = 0; p < pEnd-pStart; ++p) {
1732:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1734:     if (off >= 0) {
1735:       os = PetscMin(os, off);
1736:       oe = PetscMax(oe, off+dof);
1737:     }
1738:   }
1739:   if (start) *start = os;
1740:   if (end)   *end   = oe;
1741:   return(0);
1742: }

1744: /*@
1745:   PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields

1747:   Collective on s

1749:   Input Parameter:
1750: + s      - the PetscSection
1751: . len    - the number of subfields
1752: - fields - the subfield numbers

1754:   Output Parameter:
1755: . subs   - the subsection

1757:   Note: The section offsets now refer to a new, smaller vector.

1759:   Level: advanced

1761: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1762: @*/
1763: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1764: {
1765:   PetscInt       nF, f, c, pStart, pEnd, p, maxCdof = 0;

1769:   if (!len) return(0);
1773:   PetscSectionGetNumFields(s, &nF);
1774:   if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1775:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1776:   PetscSectionSetNumFields(*subs, len);
1777:   for (f = 0; f < len; ++f) {
1778:     const char *name   = NULL;
1779:     PetscInt   numComp = 0;

1781:     PetscSectionGetFieldName(s, fields[f], &name);
1782:     PetscSectionSetFieldName(*subs, f, name);
1783:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1784:     PetscSectionSetFieldComponents(*subs, f, numComp);
1785:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1786:       PetscSectionGetComponentName(s, fields[f], c, &name);
1787:       PetscSectionSetComponentName(*subs, f, c, name);
1788:     }
1789:   }
1790:   PetscSectionGetChart(s, &pStart, &pEnd);
1791:   PetscSectionSetChart(*subs, pStart, pEnd);
1792:   for (p = pStart; p < pEnd; ++p) {
1793:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1795:     for (f = 0; f < len; ++f) {
1796:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1797:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1798:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1799:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1800:       dof  += fdof;
1801:       cdof += cfdof;
1802:     }
1803:     PetscSectionSetDof(*subs, p, dof);
1804:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1805:     maxCdof = PetscMax(cdof, maxCdof);
1806:   }
1807:   PetscSectionSetUp(*subs);
1808:   if (maxCdof) {
1809:     PetscInt *indices;

1811:     PetscMalloc1(maxCdof, &indices);
1812:     for (p = pStart; p < pEnd; ++p) {
1813:       PetscInt cdof;

1815:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1816:       if (cdof) {
1817:         const PetscInt *oldIndices = NULL;
1818:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1820:         for (f = 0; f < len; ++f) {
1821:           PetscInt oldFoff = 0, oldf;

1823:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1824:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1825:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1826:           /* This can be sped up if we assume sorted fields */
1827:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1828:             PetscInt oldfdof = 0;
1829:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1830:             oldFoff += oldfdof;
1831:           }
1832:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1833:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1834:           numConst += cfdof;
1835:           fOff     += fdof;
1836:         }
1837:         PetscSectionSetConstraintIndices(*subs, p, indices);
1838:       }
1839:     }
1840:     PetscFree(indices);
1841:   }
1842:   return(0);
1843: }

1845: /*@
1846:   PetscSectionCreateSupersection - Create a new, larger section composed of the input sections

1848:   Collective on s

1850:   Input Parameters:
1851: + s     - the input sections
1852: - len   - the number of input sections

1854:   Output Parameter:
1855: . supers - the supersection

1857:   Note: The section offsets now refer to a new, larger vector.

1859:   Level: advanced

1861: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1862: @*/
1863: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1864: {
1865:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1869:   if (!len) return(0);
1870:   for (i = 0; i < len; ++i) {
1871:     PetscInt nf, pStarti, pEndi;

1873:     PetscSectionGetNumFields(s[i], &nf);
1874:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1875:     pStart = PetscMin(pStart, pStarti);
1876:     pEnd   = PetscMax(pEnd,   pEndi);
1877:     Nf += nf;
1878:   }
1879:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1880:   PetscSectionSetNumFields(*supers, Nf);
1881:   for (i = 0, f = 0; i < len; ++i) {
1882:     PetscInt nf, fi, ci;

1884:     PetscSectionGetNumFields(s[i], &nf);
1885:     for (fi = 0; fi < nf; ++fi, ++f) {
1886:       const char *name   = NULL;
1887:       PetscInt   numComp = 0;

1889:       PetscSectionGetFieldName(s[i], fi, &name);
1890:       PetscSectionSetFieldName(*supers, f, name);
1891:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1892:       PetscSectionSetFieldComponents(*supers, f, numComp);
1893:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1894:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1895:         PetscSectionSetComponentName(*supers, f, ci, name);
1896:       }
1897:     }
1898:   }
1899:   PetscSectionSetChart(*supers, pStart, pEnd);
1900:   for (p = pStart; p < pEnd; ++p) {
1901:     PetscInt dof = 0, cdof = 0;

1903:     for (i = 0, f = 0; i < len; ++i) {
1904:       PetscInt nf, fi, pStarti, pEndi;
1905:       PetscInt fdof = 0, cfdof = 0;

1907:       PetscSectionGetNumFields(s[i], &nf);
1908:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1909:       if ((p < pStarti) || (p >= pEndi)) continue;
1910:       for (fi = 0; fi < nf; ++fi, ++f) {
1911:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1912:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1913:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1914:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1915:         dof  += fdof;
1916:         cdof += cfdof;
1917:       }
1918:     }
1919:     PetscSectionSetDof(*supers, p, dof);
1920:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1921:     maxCdof = PetscMax(cdof, maxCdof);
1922:   }
1923:   PetscSectionSetUp(*supers);
1924:   if (maxCdof) {
1925:     PetscInt *indices;

1927:     PetscMalloc1(maxCdof, &indices);
1928:     for (p = pStart; p < pEnd; ++p) {
1929:       PetscInt cdof;

1931:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1932:       if (cdof) {
1933:         PetscInt dof, numConst = 0, fOff = 0;

1935:         for (i = 0, f = 0; i < len; ++i) {
1936:           const PetscInt *oldIndices = NULL;
1937:           PetscInt        nf, fi, pStarti, pEndi, fdof, cfdof, fc;

1939:           PetscSectionGetNumFields(s[i], &nf);
1940:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1941:           if ((p < pStarti) || (p >= pEndi)) continue;
1942:           for (fi = 0; fi < nf; ++fi, ++f) {
1943:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1944:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1945:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1946:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1947:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1948:             numConst += cfdof;
1949:           }
1950:           PetscSectionGetDof(s[i], p, &dof);
1951:           fOff += dof;
1952:         }
1953:         PetscSectionSetConstraintIndices(*supers, p, indices);
1954:       }
1955:     }
1956:     PetscFree(indices);
1957:   }
1958:   return(0);
1959: }

1961: /*@
1962:   PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh

1964:   Collective on s

1966:   Input Parameters:
1967: + s           - the PetscSection
1968: - subpointMap - a sorted list of points in the original mesh which are in the submesh

1970:   Output Parameter:
1971: . subs - the subsection

1973:   Note: The section offsets now refer to a new, smaller vector.

1975:   Level: advanced

1977: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1978: @*/
1979: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1980: {
1981:   const PetscInt *points = NULL, *indices = NULL;
1982:   PetscInt       numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;

1989:   PetscSectionGetNumFields(s, &numFields);
1990:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1991:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1992:   for (f = 0; f < numFields; ++f) {
1993:     const char *name   = NULL;
1994:     PetscInt   numComp = 0;

1996:     PetscSectionGetFieldName(s, f, &name);
1997:     PetscSectionSetFieldName(*subs, f, name);
1998:     PetscSectionGetFieldComponents(s, f, &numComp);
1999:     PetscSectionSetFieldComponents(*subs, f, numComp);
2000:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2001:       PetscSectionGetComponentName(s, f, c, &name);
2002:       PetscSectionSetComponentName(*subs, f, c, name);
2003:     }
2004:   }
2005:   /* For right now, we do not try to squeeze the subchart */
2006:   if (subpointMap) {
2007:     ISGetSize(subpointMap, &numSubpoints);
2008:     ISGetIndices(subpointMap, &points);
2009:   }
2010:   PetscSectionGetChart(s, &pStart, &pEnd);
2011:   PetscSectionSetChart(*subs, 0, numSubpoints);
2012:   for (p = pStart; p < pEnd; ++p) {
2013:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

2015:     PetscFindInt(p, numSubpoints, points, &subp);
2016:     if (subp < 0) continue;
2017:     for (f = 0; f < numFields; ++f) {
2018:       PetscSectionGetFieldDof(s, p, f, &fdof);
2019:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
2020:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
2021:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
2022:     }
2023:     PetscSectionGetDof(s, p, &dof);
2024:     PetscSectionSetDof(*subs, subp, dof);
2025:     PetscSectionGetConstraintDof(s, p, &cdof);
2026:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
2027:   }
2028:   PetscSectionSetUp(*subs);
2029:   /* Change offsets to original offsets */
2030:   for (p = pStart; p < pEnd; ++p) {
2031:     PetscInt off, foff = 0;

2033:     PetscFindInt(p, numSubpoints, points, &subp);
2034:     if (subp < 0) continue;
2035:     for (f = 0; f < numFields; ++f) {
2036:       PetscSectionGetFieldOffset(s, p, f, &foff);
2037:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
2038:     }
2039:     PetscSectionGetOffset(s, p, &off);
2040:     PetscSectionSetOffset(*subs, subp, off);
2041:   }
2042:   /* Copy constraint indices */
2043:   for (subp = 0; subp < numSubpoints; ++subp) {
2044:     PetscInt cdof;

2046:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
2047:     if (cdof) {
2048:       for (f = 0; f < numFields; ++f) {
2049:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
2050:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2051:       }
2052:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
2053:       PetscSectionSetConstraintIndices(*subs, subp, indices);
2054:     }
2055:   }
2056:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2057:   return(0);
2058: }

2060: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2061: {
2062:   PetscInt       p;
2063:   PetscMPIInt    rank;

2067:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2068:   PetscViewerASCIIPushSynchronized(viewer);
2069:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2070:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2071:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2072:       PetscInt b;

2074:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2075:       if (s->bcIndices) {
2076:         for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2077:           PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2078:         }
2079:       }
2080:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2081:     } else {
2082:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2083:     }
2084:   }
2085:   PetscViewerFlush(viewer);
2086:   PetscViewerASCIIPopSynchronized(viewer);
2087:   if (s->sym) {
2088:     PetscViewerASCIIPushTab(viewer);
2089:     PetscSectionSymView(s->sym,viewer);
2090:     PetscViewerASCIIPopTab(viewer);
2091:   }
2092:   return(0);
2093: }

2095: /*@C
2096:    PetscSectionViewFromOptions - View from Options

2098:    Collective on PetscSection

2100:    Input Parameters:
2101: +  A - the PetscSection object to view
2102: .  obj - Optional object
2103: -  name - command line option

2105:    Level: intermediate
2106: .seealso:  PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2107: @*/
2108: PetscErrorCode  PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2109: {

2114:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
2115:   return(0);
2116: }

2118: /*@C
2119:   PetscSectionView - Views a PetscSection

2121:   Collective on PetscSection

2123:   Input Parameters:
2124: + s - the PetscSection object to view
2125: - v - the viewer

2127:   Level: beginner

2129: .seealso PetscSectionCreate(), PetscSectionDestroy()
2130: @*/
2131: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2132: {
2133:   PetscBool      isascii;
2134:   PetscInt       f;

2139:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2141:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2142:   if (isascii) {
2143:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2144:     if (s->numFields) {
2145:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2146:       for (f = 0; f < s->numFields; ++f) {
2147:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
2148:         PetscSectionView_ASCII(s->field[f], viewer);
2149:       }
2150:     } else {
2151:       PetscSectionView_ASCII(s, viewer);
2152:     }
2153:   }
2154:   return(0);
2155: }

2157: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2158: {
2160:   PetscSectionClosurePermVal clVal;

2163:   if (!section->clHash) return(0);
2164:   kh_foreach_value(section->clHash, clVal, {
2165:       PetscFree(clVal.perm);
2166:       PetscFree(clVal.invPerm);
2167:     });
2168:   kh_destroy(ClPerm, section->clHash);
2169:   section->clHash = NULL;
2170:   return(0);
2171: }

2173: /*@
2174:   PetscSectionReset - Frees all section data.

2176:   Not collective

2178:   Input Parameters:
2179: . s - the PetscSection

2181:   Level: beginner

2183: .seealso: PetscSection, PetscSectionCreate()
2184: @*/
2185: PetscErrorCode PetscSectionReset(PetscSection s)
2186: {
2187:   PetscInt       f, c;

2192:   for (f = 0; f < s->numFields; ++f) {
2193:     PetscSectionDestroy(&s->field[f]);
2194:     PetscFree(s->fieldNames[f]);
2195:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2196:       PetscFree(s->compNames[f][c]);
2197:     }
2198:     PetscFree(s->compNames[f]);
2199:   }
2200:   PetscFree(s->numFieldComponents);
2201:   PetscFree(s->fieldNames);
2202:   PetscFree(s->compNames);
2203:   PetscFree(s->field);
2204:   PetscSectionDestroy(&s->bc);
2205:   PetscFree(s->bcIndices);
2206:   PetscFree2(s->atlasDof, s->atlasOff);
2207:   PetscSectionDestroy(&s->clSection);
2208:   ISDestroy(&s->clPoints);
2209:   ISDestroy(&s->perm);
2210:   PetscSectionResetClosurePermutation(s);
2211:   PetscSectionSymDestroy(&s->sym);
2212:   PetscSectionDestroy(&s->clSection);
2213:   ISDestroy(&s->clPoints);

2215:   s->pStart    = -1;
2216:   s->pEnd      = -1;
2217:   s->maxDof    = 0;
2218:   s->setup     = PETSC_FALSE;
2219:   s->numFields = 0;
2220:   s->clObj     = NULL;
2221:   return(0);
2222: }

2224: /*@
2225:   PetscSectionDestroy - Frees a section object and frees its range if that exists.

2227:   Not collective

2229:   Input Parameters:
2230: . s - the PetscSection

2232:   Level: beginner

2234: .seealso: PetscSection, PetscSectionCreate()
2235: @*/
2236: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2237: {

2241:   if (!*s) return(0);
2243:   if (--((PetscObject)(*s))->refct > 0) {
2244:     *s = NULL;
2245:     return(0);
2246:   }
2247:   PetscSectionReset(*s);
2248:   PetscHeaderDestroy(s);
2249:   return(0);
2250: }

2252: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2253: {
2254:   const PetscInt p = point - s->pStart;

2258:   *values = &baseArray[s->atlasOff[p]];
2259:   return(0);
2260: }

2262: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2263: {
2264:   PetscInt       *array;
2265:   const PetscInt p           = point - s->pStart;
2266:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2267:   PetscInt       cDim        = 0;

2272:   PetscSectionGetConstraintDof(s, p, &cDim);
2273:   array = &baseArray[s->atlasOff[p]];
2274:   if (!cDim) {
2275:     if (orientation >= 0) {
2276:       const PetscInt dim = s->atlasDof[p];
2277:       PetscInt       i;

2279:       if (mode == INSERT_VALUES) {
2280:         for (i = 0; i < dim; ++i) array[i] = values[i];
2281:       } else {
2282:         for (i = 0; i < dim; ++i) array[i] += values[i];
2283:       }
2284:     } else {
2285:       PetscInt offset = 0;
2286:       PetscInt j      = -1, field, i;

2288:       for (field = 0; field < s->numFields; ++field) {
2289:         const PetscInt dim = s->field[field]->atlasDof[p];

2291:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2292:         offset += dim;
2293:       }
2294:     }
2295:   } else {
2296:     if (orientation >= 0) {
2297:       const PetscInt dim  = s->atlasDof[p];
2298:       PetscInt       cInd = 0, i;
2299:       const PetscInt *cDof;

2301:       PetscSectionGetConstraintIndices(s, point, &cDof);
2302:       if (mode == INSERT_VALUES) {
2303:         for (i = 0; i < dim; ++i) {
2304:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2305:           array[i] = values[i];
2306:         }
2307:       } else {
2308:         for (i = 0; i < dim; ++i) {
2309:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2310:           array[i] += values[i];
2311:         }
2312:       }
2313:     } else {
2314:       const PetscInt *cDof;
2315:       PetscInt       offset  = 0;
2316:       PetscInt       cOffset = 0;
2317:       PetscInt       j       = 0, field;

2319:       PetscSectionGetConstraintIndices(s, point, &cDof);
2320:       for (field = 0; field < s->numFields; ++field) {
2321:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2322:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2323:         const PetscInt sDim = dim - tDim;
2324:         PetscInt       cInd = 0, i,k;

2326:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2327:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2328:           array[j] = values[k];
2329:         }
2330:         offset  += dim;
2331:         cOffset += dim - tDim;
2332:       }
2333:     }
2334:   }
2335:   return(0);
2336: }

2338: /*@C
2339:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2341:   Not collective

2343:   Input Parameter:
2344: . s - The PetscSection

2346:   Output Parameter:
2347: . hasConstraints - flag indicating that the section has constrained dofs

2349:   Level: intermediate

2351: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2352: @*/
2353: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2354: {
2358:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2359:   return(0);
2360: }

2362: /*@C
2363:   PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained

2365:   Not collective

2367:   Input Parameters:
2368: + s     - The PetscSection
2369: - point - The point

2371:   Output Parameter:
2372: . indices - The constrained dofs

2374:   Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()

2376:   Level: intermediate

2378: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2379: @*/
2380: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2381: {

2386:   if (s->bc) {
2387:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2388:   } else *indices = NULL;
2389:   return(0);
2390: }

2392: /*@C
2393:   PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained

2395:   Not collective

2397:   Input Parameters:
2398: + s     - The PetscSection
2399: . point - The point
2400: - indices - The constrained dofs

2402:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2404:   Level: intermediate

2406: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2407: @*/
2408: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2409: {

2414:   if (s->bc) {
2415:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2416:   }
2417:   return(0);
2418: }

2420: /*@C
2421:   PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained

2423:   Not collective

2425:   Input Parameters:
2426: + s     - The PetscSection
2427: . field  - The field number
2428: - point - The point

2430:   Output Parameter:
2431: . indices - The constrained dofs sorted in ascending order

2433:   Notes:
2434:   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().

2436:   Fortran Note:
2437:   In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()

2439:   Level: intermediate

2441: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2442: @*/
2443: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2444: {

2449:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
2450:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2451:   return(0);
2452: }

2454: /*@C
2455:   PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained

2457:   Not collective

2459:   Input Parameters:
2460: + s       - The PetscSection
2461: . point   - The point
2462: . field   - The field number
2463: - indices - The constrained dofs

2465:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2467:   Level: intermediate

2469: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2470: @*/
2471: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2472: {

2477:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
2478:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2479:   return(0);
2480: }

2482: /*@
2483:   PetscSectionPermute - Reorder the section according to the input point permutation

2485:   Collective on PetscSection

2487:   Input Parameter:
2488: + section - The PetscSection object
2489: - perm - The point permutation, old point p becomes new point perm[p]

2491:   Output Parameter:
2492: . sectionNew - The permuted PetscSection

2494:   Level: intermediate

2496: .seealso: MatPermute()
2497: @*/
2498: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2499: {
2500:   PetscSection    s = section, sNew;
2501:   const PetscInt *perm;
2502:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;
2503:   PetscErrorCode  ierr;

2509:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2510:   PetscSectionGetNumFields(s, &numFields);
2511:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2512:   for (f = 0; f < numFields; ++f) {
2513:     const char *name;
2514:     PetscInt    numComp;

2516:     PetscSectionGetFieldName(s, f, &name);
2517:     PetscSectionSetFieldName(sNew, f, name);
2518:     PetscSectionGetFieldComponents(s, f, &numComp);
2519:     PetscSectionSetFieldComponents(sNew, f, numComp);
2520:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2521:       PetscSectionGetComponentName(s, f, c, &name);
2522:       PetscSectionSetComponentName(sNew, f, c, name);
2523:     }
2524:   }
2525:   ISGetLocalSize(permutation, &numPoints);
2526:   ISGetIndices(permutation, &perm);
2527:   PetscSectionGetChart(s, &pStart, &pEnd);
2528:   PetscSectionSetChart(sNew, pStart, pEnd);
2529:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2530:   for (p = pStart; p < pEnd; ++p) {
2531:     PetscInt dof, cdof;

2533:     PetscSectionGetDof(s, p, &dof);
2534:     PetscSectionSetDof(sNew, perm[p], dof);
2535:     PetscSectionGetConstraintDof(s, p, &cdof);
2536:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2537:     for (f = 0; f < numFields; ++f) {
2538:       PetscSectionGetFieldDof(s, p, f, &dof);
2539:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2540:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2541:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2542:     }
2543:   }
2544:   PetscSectionSetUp(sNew);
2545:   for (p = pStart; p < pEnd; ++p) {
2546:     const PetscInt *cind;
2547:     PetscInt        cdof;

2549:     PetscSectionGetConstraintDof(s, p, &cdof);
2550:     if (cdof) {
2551:       PetscSectionGetConstraintIndices(s, p, &cind);
2552:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2553:     }
2554:     for (f = 0; f < numFields; ++f) {
2555:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2556:       if (cdof) {
2557:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2558:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2559:       }
2560:     }
2561:   }
2562:   ISRestoreIndices(permutation, &perm);
2563:   *sectionNew = sNew;
2564:   return(0);
2565: }

2567: /*@
2568:   PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section

2570:   Collective on section

2572:   Input Parameters:
2573: + section   - The PetscSection
2574: . obj       - A PetscObject which serves as the key for this index
2575: . clSection - Section giving the size of the closure of each point
2576: - clPoints  - IS giving the points in each closure

2578:   Note: We compress out closure points with no dofs in this section

2580:   Level: advanced

2582: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2583: @*/
2584: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2585: {

2592:   if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
2593:   section->clObj     = obj;
2594:   PetscObjectReference((PetscObject)clSection);
2595:   PetscObjectReference((PetscObject)clPoints);
2596:   PetscSectionDestroy(&section->clSection);
2597:   ISDestroy(&section->clPoints);
2598:   section->clSection = clSection;
2599:   section->clPoints  = clPoints;
2600:   return(0);
2601: }

2603: /*@
2604:   PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section

2606:   Collective on section

2608:   Input Parameters:
2609: + section   - The PetscSection
2610: - obj       - A PetscObject which serves as the key for this index

2612:   Output Parameters:
2613: + clSection - Section giving the size of the closure of each point
2614: - clPoints  - IS giving the points in each closure

2616:   Note: We compress out closure points with no dofs in this section

2618:   Level: advanced

2620: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2621: @*/
2622: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2623: {
2625:   if (section->clObj == obj) {
2626:     if (clSection) *clSection = section->clSection;
2627:     if (clPoints)  *clPoints  = section->clPoints;
2628:   } else {
2629:     if (clSection) *clSection = NULL;
2630:     if (clPoints)  *clPoints  = NULL;
2631:   }
2632:   return(0);
2633: }

2635: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2636: {
2637:   PetscInt       i;
2638:   khiter_t iter;
2639:   int new_entry;
2640:   PetscSectionClosurePermKey key = {depth, clSize};
2641:   PetscSectionClosurePermVal *val;

2645:   if (section->clObj != obj) {
2646:     PetscSectionDestroy(&section->clSection);
2647:     ISDestroy(&section->clPoints);
2648:   }
2649:   section->clObj = obj;
2650:   if (!section->clHash) {PetscClPermCreate(&section->clHash);}
2651:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2652:   val = &kh_val(section->clHash, iter);
2653:   if (!new_entry) {
2654:     PetscFree(val->perm);
2655:     PetscFree(val->invPerm);
2656:   }
2657:   if (mode == PETSC_COPY_VALUES) {
2658:     PetscMalloc1(clSize, &val->perm);
2659:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2660:     PetscArraycpy(val->perm, clPerm, clSize);
2661:   } else if (mode == PETSC_OWN_POINTER) {
2662:     val->perm = clPerm;
2663:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2664:   PetscMalloc1(clSize, &val->invPerm);
2665:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2666:   return(0);
2667: }

2669: /*@
2670:   PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2672:   Not Collective

2674:   Input Parameters:
2675: + section - The PetscSection
2676: . obj     - A PetscObject which serves as the key for this index (usually a DM)
2677: . depth   - Depth of points on which to apply the given permutation
2678: - perm    - Permutation of the cell dof closure

2680:   Note:
2681:   The specified permutation will only be applied to points at depth whose closure size matches the length of perm.  In a
2682:   mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2683:   each topology and degree.

2685:   This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2686:   exotic/enriched spaces on mixed topology meshes.

2688:   Level: intermediate

2690: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2691: @*/
2692: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2693: {
2694:   const PetscInt *clPerm = NULL;
2695:   PetscInt        clSize = 0;
2696:   PetscErrorCode  ierr;

2699:   if (perm) {
2700:     ISGetLocalSize(perm, &clSize);
2701:     ISGetIndices(perm, &clPerm);
2702:   }
2703:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2704:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2705:   return(0);
2706: }

2708: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2709: {

2713:   if (section->clObj == obj) {
2714:     PetscSectionClosurePermKey k = {depth, size};
2715:     PetscSectionClosurePermVal v;
2716:     PetscClPermGet(section->clHash, k, &v);
2717:     if (perm) *perm = v.perm;
2718:   } else {
2719:     if (perm) *perm = NULL;
2720:   }
2721:   return(0);
2722: }

2724: /*@
2725:   PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2727:   Not collective

2729:   Input Parameters:
2730: + section   - The PetscSection
2731: . obj       - A PetscObject which serves as the key for this index (usually a DM)
2732: . depth     - Depth stratum on which to obtain closure permutation
2733: - clSize    - Closure size to be permuted (e.g., may vary with element topology and degree)

2735:   Output Parameter:
2736: . perm - The dof closure permutation

2738:   Note:
2739:   The user must destroy the IS that is returned.

2741:   Level: intermediate

2743: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2744: @*/
2745: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2746: {
2747:   const PetscInt *clPerm;
2748:   PetscErrorCode  ierr;

2751:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2752:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2753:   return(0);
2754: }

2756: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2757: {

2761:   if (section->clObj == obj && section->clHash) {
2762:     PetscSectionClosurePermKey k = {depth, size};
2763:     PetscSectionClosurePermVal v;
2764:     PetscClPermGet(section->clHash, k, &v);
2765:     if (perm) *perm = v.invPerm;
2766:   } else {
2767:     if (perm) *perm = NULL;
2768:   }
2769:   return(0);
2770: }

2772: /*@
2773:   PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.

2775:   Not collective

2777:   Input Parameters:
2778: + section   - The PetscSection
2779: . obj       - A PetscObject which serves as the key for this index (usually a DM)
2780: . depth     - Depth stratum on which to obtain closure permutation
2781: - clSize    - Closure size to be permuted (e.g., may vary with element topology and degree)

2783:   Output Parameters:
2784: . perm - The dof closure permutation

2786:   Note:
2787:   The user must destroy the IS that is returned.

2789:   Level: intermediate

2791: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2792: @*/
2793: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2794: {
2795:   const PetscInt *clPerm;
2796:   PetscErrorCode  ierr;

2799:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2800:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2801:   return(0);
2802: }

2804: /*@
2805:   PetscSectionGetField - Get the subsection associated with a single field

2807:   Input Parameters:
2808: + s     - The PetscSection
2809: - field - The field number

2811:   Output Parameter:
2812: . subs  - The subsection for the given field

2814:   Level: intermediate

2816: .seealso: PetscSectionSetNumFields()
2817: @*/
2818: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2819: {
2823:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
2824:   *subs = s->field[field];
2825:   return(0);
2826: }

2828: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2829: PetscFunctionList PetscSectionSymList = NULL;

2831: /*@
2832:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2834:   Collective

2836:   Input Parameter:
2837: . comm - the MPI communicator

2839:   Output Parameter:
2840: . sym - pointer to the new set of symmetries

2842:   Level: developer

2844: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2845: @*/
2846: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2847: {

2852:   ISInitializePackage();
2853:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2854:   return(0);
2855: }

2857: /*@C
2858:   PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.

2860:   Collective on PetscSectionSym

2862:   Input Parameters:
2863: + sym    - The section symmetry object
2864: - method - The name of the section symmetry type

2866:   Level: developer

2868: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2869: @*/
2870: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2871: {
2872:   PetscErrorCode (*r)(PetscSectionSym);
2873:   PetscBool      match;

2878:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2879:   if (match) return(0);

2881:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2882:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2883:   if (sym->ops->destroy) {
2884:     (*sym->ops->destroy)(sym);
2885:     sym->ops->destroy = NULL;
2886:   }
2887:   (*r)(sym);
2888:   PetscObjectChangeTypeName((PetscObject)sym,method);
2889:   return(0);
2890: }


2893: /*@C
2894:   PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.

2896:   Not Collective

2898:   Input Parameter:
2899: . sym  - The section symmetry

2901:   Output Parameter:
2902: . type - The index set type name

2904:   Level: developer

2906: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2907: @*/
2908: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2909: {
2913:   *type = ((PetscObject)sym)->type_name;
2914:   return(0);
2915: }

2917: /*@C
2918:   PetscSectionSymRegister - Adds a new section symmetry implementation

2920:   Not Collective

2922:   Input Parameters:
2923: + name        - The name of a new user-defined creation routine
2924: - create_func - The creation routine itself

2926:   Notes:
2927:   PetscSectionSymRegister() may be called multiple times to add several user-defined vectors

2929:   Level: developer

2931: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2932: @*/
2933: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2934: {

2938:   ISInitializePackage();
2939:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2940:   return(0);
2941: }

2943: /*@
2944:    PetscSectionSymDestroy - Destroys a section symmetry.

2946:    Collective on PetscSectionSym

2948:    Input Parameters:
2949: .  sym - the section symmetry

2951:    Level: developer

2953: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2954: @*/
2955: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2956: {
2957:   SymWorkLink    link,next;

2961:   if (!*sym) return(0);
2963:   if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
2964:   if ((*sym)->ops->destroy) {
2965:     (*(*sym)->ops->destroy)(*sym);
2966:   }
2967:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2968:   for (link=(*sym)->workin; link; link=next) {
2969:     next = link->next;
2970:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2971:     PetscFree(link);
2972:   }
2973:   (*sym)->workin = NULL;
2974:   PetscHeaderDestroy(sym);
2975:   return(0);
2976: }

2978: /*@C
2979:    PetscSectionSymView - Displays a section symmetry

2981:    Collective on PetscSectionSym

2983:    Input Parameters:
2984: +  sym - the index set
2985: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2987:    Level: developer

2989: .seealso: PetscViewerASCIIOpen()
2990: @*/
2991: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2992: {

2997:   if (!viewer) {
2998:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2999:   }
3002:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3003:   if (sym->ops->view) {
3004:     (*sym->ops->view)(sym,viewer);
3005:   }
3006:   return(0);
3007: }

3009: /*@
3010:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3012:   Collective on PetscSection

3014:   Input Parameters:
3015: + section - the section describing data layout
3016: - sym - the symmetry describing the affect of orientation on the access of the data

3018:   Level: developer

3020: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3021: @*/
3022: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3023: {

3028:   PetscSectionSymDestroy(&(section->sym));
3029:   if (sym) {
3032:     PetscObjectReference((PetscObject) sym);
3033:   }
3034:   section->sym = sym;
3035:   return(0);
3036: }

3038: /*@
3039:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3041:   Not collective

3043:   Input Parameters:
3044: . section - the section describing data layout

3046:   Output Parameters:
3047: . sym - the symmetry describing the affect of orientation on the access of the data

3049:   Level: developer

3051: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3052: @*/
3053: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3054: {
3057:   *sym = section->sym;
3058:   return(0);
3059: }

3061: /*@
3062:   PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section

3064:   Collective on PetscSection

3066:   Input Parameters:
3067: + section - the section describing data layout
3068: . field - the field number
3069: - sym - the symmetry describing the affect of orientation on the access of the data

3071:   Level: developer

3073: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3074: @*/
3075: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3076: {

3081:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
3082:   PetscSectionSetSym(section->field[field],sym);
3083:   return(0);
3084: }

3086: /*@
3087:   PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section

3089:   Collective on PetscSection

3091:   Input Parameters:
3092: + section - the section describing data layout
3093: - field - the field number

3095:   Output Parameters:
3096: . sym - the symmetry describing the affect of orientation on the access of the data

3098:   Level: developer

3100: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3101: @*/
3102: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3103: {
3106:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
3107:   *sym = section->field[field]->sym;
3108:   return(0);
3109: }

3111: /*@C
3112:   PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.

3114:   Not collective

3116:   Input Parameters:
3117: + section - the section
3118: . numPoints - the number of points
3119: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3120:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3121:     context, see DMPlexGetConeOrientation()).

3123:   Output Parameter:
3124: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3125: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3126:     identity).

3128:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3129: .vb
3130:      const PetscInt    **perms;
3131:      const PetscScalar **rots;
3132:      PetscInt            lOffset;

3134:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3135:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3136:        PetscInt           point = points[2*i], dof, sOffset;
3137:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3138:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3140:        PetscSectionGetDof(section,point,&dof);
3141:        PetscSectionGetOffset(section,point,&sOffset);

3143:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3144:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3145:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3146:        lOffset += dof;
3147:      }
3148:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3149: .ve

3151:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3152: .vb
3153:      const PetscInt    **perms;
3154:      const PetscScalar **rots;
3155:      PetscInt            lOffset;

3157:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3158:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3159:        PetscInt           point = points[2*i], dof, sOffset;
3160:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3161:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3163:        PetscSectionGetDof(section,point,&dof);
3164:        PetscSectionGetOffset(section,point,&sOff);

3166:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3167:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3168:        offset += dof;
3169:      }
3170:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3171: .ve

3173:   Level: developer

3175: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3176: @*/
3177: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3178: {
3179:   PetscSectionSym sym;
3180:   PetscErrorCode  ierr;

3185:   if (perms) *perms = NULL;
3186:   if (rots)  *rots  = NULL;
3187:   sym = section->sym;
3188:   if (sym && (perms || rots)) {
3189:     SymWorkLink link;

3191:     if (sym->workin) {
3192:       link        = sym->workin;
3193:       sym->workin = sym->workin->next;
3194:     } else {
3195:       PetscNewLog(sym,&link);
3196:     }
3197:     if (numPoints > link->numPoints) {
3198:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3199:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3200:       link->numPoints = numPoints;
3201:     }
3202:     link->next   = sym->workout;
3203:     sym->workout = link;
3204:     PetscArrayzero((PetscInt**)link->perms,numPoints);
3205:     PetscArrayzero((PetscInt**)link->rots,numPoints);
3206:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3207:     if (perms) *perms = link->perms;
3208:     if (rots)  *rots  = link->rots;
3209:   }
3210:   return(0);
3211: }

3213: /*@C
3214:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3216:   Not collective

3218:   Input Parameters:
3219: + section - the section
3220: . numPoints - the number of points
3221: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3222:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3223:     context, see DMPlexGetConeOrientation()).

3225:   Output Parameter:
3226: + perms - The permutations for the given orientations: set to NULL at conclusion
3227: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3229:   Level: developer

3231: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3232: @*/
3233: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3234: {
3235:   PetscSectionSym sym;

3239:   sym = section->sym;
3240:   if (sym && (perms || rots)) {
3241:     SymWorkLink *p,link;

3243:     for (p=&sym->workout; (link=*p); p=&link->next) {
3244:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3245:         *p          = link->next;
3246:         link->next  = sym->workin;
3247:         sym->workin = link;
3248:         if (perms) *perms = NULL;
3249:         if (rots)  *rots  = NULL;
3250:         return(0);
3251:       }
3252:     }
3253:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3254:   }
3255:   return(0);
3256: }

3258: /*@C
3259:   PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.

3261:   Not collective

3263:   Input Parameters:
3264: + section - the section
3265: . field - the field of the section
3266: . numPoints - the number of points
3267: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3268:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3269:     context, see DMPlexGetConeOrientation()).

3271:   Output Parameter:
3272: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3273: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3274:     identity).

3276:   Level: developer

3278: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3279: @*/
3280: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3281: {

3286:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
3287:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3288:   return(0);
3289: }

3291: /*@C
3292:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3294:   Not collective

3296:   Input Parameters:
3297: + section - the section
3298: . field - the field number
3299: . numPoints - the number of points
3300: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3301:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3302:     context, see DMPlexGetConeOrientation()).

3304:   Output Parameter:
3305: + perms - The permutations for the given orientations: set to NULL at conclusion
3306: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3308:   Level: developer

3310: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3311: @*/
3312: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3313: {

3318:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
3319:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3320:   return(0);
3321: }

3323: /*@
3324:   PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset

3326:   Not collective

3328:   Input Parameter:
3329: . s - the global PetscSection

3331:   Output Parameters:
3332: . flg - the flag

3334:   Level: developer

3336: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3337: @*/
3338: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3339: {
3342:   *flg = s->useFieldOff;
3343:   return(0);
3344: }

3346: /*@
3347:   PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset

3349:   Not collective

3351:   Input Parameters:
3352: + s   - the global PetscSection
3353: - flg - the flag

3355:   Level: developer

3357: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3358: @*/
3359: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3360: {
3363:   s->useFieldOff = flg;
3364:   return(0);
3365: }

3367: #define PetscSectionExpandPoints_Loop(TYPE) \
3368: { \
3369:   PetscInt i, n, o0, o1, size; \
3370:   TYPE *a0 = (TYPE*)origArray, *a1; \
3371:   PetscSectionGetStorageSize(s, &size); \
3372:   PetscMalloc1(size, &a1); \
3373:   for (i=0; i<npoints; i++) { \
3374:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3375:     PetscSectionGetOffset(s, i, &o1); \
3376:     PetscSectionGetDof(s, i, &n); \
3377:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3378:   } \
3379:   *newArray = (void*)a1; \
3380: }

3382: /*@
3383:   PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.

3385:   Not collective

3387:   Input Parameters:
3388: + origSection - the PetscSection describing the layout of the array
3389: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3390: . origArray - the array; its size must be equal to the storage size of origSection
3391: - points - IS with points to extract; its indices must lie in the chart of origSection

3393:   Output Parameters:
3394: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3395: - newArray - the array of the extracted DOFs; its size is the storage size of newSection

3397:   Level: developer

3399: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3400: @*/
3401: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3402: {
3403:   PetscSection        s;
3404:   const PetscInt      *points_;
3405:   PetscInt            i, n, npoints, pStart, pEnd;
3406:   PetscMPIInt         unitsize;
3407:   PetscErrorCode      ierr;

3415:   MPI_Type_size(dataType, &unitsize);
3416:   ISGetLocalSize(points, &npoints);
3417:   ISGetIndices(points, &points_);
3418:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3419:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3420:   PetscSectionSetChart(s, 0, npoints);
3421:   for (i=0; i<npoints; i++) {
3422:     if (PetscUnlikely(points_[i] < pStart || points_[i] >= pEnd)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %d (index %d) in input IS out of input section's chart", points_[i], i);
3423:     PetscSectionGetDof(origSection, points_[i], &n);
3424:     PetscSectionSetDof(s, i, n);
3425:   }
3426:   PetscSectionSetUp(s);
3427:   if (newArray) {
3428:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3429:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3430:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3431:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3432:   }
3433:   if (newSection) {
3434:     *newSection = s;
3435:   } else {
3436:     PetscSectionDestroy(&s);
3437:   }
3438:   return(0);
3439: }