Actual source code: section.c

petsc-master 2020-01-17
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)->maxDof             = 0;
 53:   (*s)->atlasDof           = NULL;
 54:   (*s)->atlasOff           = NULL;
 55:   (*s)->bc                 = NULL;
 56:   (*s)->bcIndices          = NULL;
 57:   (*s)->setup              = PETSC_FALSE;
 58:   (*s)->numFields          = 0;
 59:   (*s)->fieldNames         = NULL;
 60:   (*s)->field              = NULL;
 61:   (*s)->useFieldOff        = PETSC_FALSE;
 62:   (*s)->clObj              = NULL;
 63:   (*s)->clSection          = NULL;
 64:   (*s)->clPoints           = NULL;
 65:   (*s)->clSize             = 0;
 66:   (*s)->clPerm             = NULL;
 67:   (*s)->clInvPerm          = 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, 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 *name   = NULL;
101:     PetscInt   numComp = 0;

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

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

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

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

159:   Collective

161:   Input Parameter:
162: . section - the PetscSection

164:   Output Parameter:
165: . newSection - the copy

167:   Level: beginner

169: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
170: @*/
171: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
172: {

178:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
179:   PetscSectionCopy(section, *newSection);
180:   return(0);
181: }

183: /*@
184:   PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database

186:   Collective on PetscSection

188:   Input Parameter:
189: . section - the PetscSection

191:   Options Database:
192: . -petscsection_point_major the dof order

194:   Level: intermediate

196: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
197: @*/
198: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
199: {

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

213: /*@
214:   PetscSectionCompare - Compares two sections

216:   Collective on PetscSection

218:   Input Parameterss:
219: + s1 - the first PetscSection
220: - s2 - the second PetscSection

222:   Output Parameter:
223: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise

225:   Level: intermediate

227:   Notes:
228:   Field names are disregarded.

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

245:   flg = PETSC_FALSE;

247:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
248:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
249:     *congruent = PETSC_FALSE;
250:     return(0);
251:   }

253:   PetscSectionGetChart(s1, &pStart, &pEnd);
254:   PetscSectionGetChart(s2, &n1, &n2);
255:   if (pStart != n1 || pEnd != n2) goto not_congruent;

257:   PetscSectionGetPermutation(s1, &perm1);
258:   PetscSectionGetPermutation(s2, &perm2);
259:   if (perm1 && perm2) {
260:     ISEqual(perm1, perm2, congruent);
261:     if (!(*congruent)) goto not_congruent;
262:   } else if (perm1 != perm2) goto not_congruent;

264:   for (p = pStart; p < pEnd; ++p) {
265:     PetscSectionGetOffset(s1, p, &n1);
266:     PetscSectionGetOffset(s2, p, &n2);
267:     if (n1 != n2) goto not_congruent;

269:     PetscSectionGetDof(s1, p, &n1);
270:     PetscSectionGetDof(s2, p, &n2);
271:     if (n1 != n2) goto not_congruent;

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

277:     PetscSectionGetConstraintIndices(s1, p, &idx1);
278:     PetscSectionGetConstraintIndices(s2, p, &idx2);
279:     PetscArraycmp(idx1, idx2, ncdof, congruent);
280:     if (!(*congruent)) goto not_congruent;
281:   }

283:   PetscSectionGetNumFields(s1, &nfields);
284:   PetscSectionGetNumFields(s2, &n2);
285:   if (nfields != n2) goto not_congruent;

287:   for (f = 0; f < nfields; ++f) {
288:     PetscSectionGetFieldComponents(s1, f, &n1);
289:     PetscSectionGetFieldComponents(s2, f, &n2);
290:     if (n1 != n2) goto not_congruent;

292:     for (p = pStart; p < pEnd; ++p) {
293:       PetscSectionGetFieldOffset(s1, p, f, &n1);
294:       PetscSectionGetFieldOffset(s2, p, f, &n2);
295:       if (n1 != n2) goto not_congruent;

297:       PetscSectionGetFieldDof(s1, p, f, &n1);
298:       PetscSectionGetFieldDof(s2, p, f, &n2);
299:       if (n1 != n2) goto not_congruent;

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

305:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
306:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
307:       PetscArraycmp(idx1, idx2, nfcdof, congruent);
308:       if (!(*congruent)) goto not_congruent;
309:     }
310:   }

312:   flg = PETSC_TRUE;
313: not_congruent:
314:   MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
315:   return(0);
316: }

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

321:   Not collective

323:   Input Parameter:
324: . s - the PetscSection

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

329:   Level: intermediate

331: .seealso: PetscSectionSetNumFields()
332: @*/
333: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
334: {
338:   *numFields = s->numFields;
339:   return(0);
340: }

342: /*@
343:   PetscSectionSetNumFields - Sets the number of fields.

345:   Not collective

347:   Input Parameters:
348: + s - the PetscSection
349: - numFields - the number of fields

351:   Level: intermediate

353: .seealso: PetscSectionGetNumFields()
354: @*/
355: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
356: {
357:   PetscInt       f;

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

365:   s->numFields = numFields;
366:   PetscMalloc1(s->numFields, &s->numFieldComponents);
367:   PetscMalloc1(s->numFields, &s->fieldNames);
368:   PetscMalloc1(s->numFields, &s->field);
369:   for (f = 0; f < s->numFields; ++f) {
370:     char name[64];

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

374:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
375:     PetscSNPrintf(name, 64, "Field_%D", f);
376:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
377:   }
378:   return(0);
379: }

381: /*@C
382:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

384:   Not Collective

386:   Input Parameters:
387: + s     - the PetscSection
388: - field - the field number

390:   Output Parameter:
391: . fieldName - the field name

393:   Level: intermediate

395: .seealso: PetscSectionSetFieldName()
396: @*/
397: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
398: {
402:   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);
403:   *fieldName = s->fieldNames[field];
404:   return(0);
405: }

407: /*@C
408:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

410:   Not Collective

412:   Input Parameters:
413: + s     - the PetscSection
414: . field - the field number
415: - fieldName - the field name

417:   Level: intermediate

419: .seealso: PetscSectionGetFieldName()
420: @*/
421: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
422: {

428:   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);
429:   PetscFree(s->fieldNames[field]);
430:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
431:   return(0);
432: }

434: /*@
435:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

437:   Not collective

439:   Input Parameters:
440: + s - the PetscSection
441: - field - the field number

443:   Output Parameter:
444: . numComp - the number of field components

446:   Level: intermediate

448: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
449: @*/
450: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
451: {
455:   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);
456:   *numComp = s->numFieldComponents[field];
457:   return(0);
458: }

460: /*@
461:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

463:   Not collective

465:   Input Parameters:
466: + s - the PetscSection
467: . field - the field number
468: - numComp - the number of field components

470:   Level: intermediate

472: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
473: @*/
474: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
475: {
478:   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);
479:   s->numFieldComponents[field] = numComp;
480:   return(0);
481: }

483: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
484: {

488:   if (!s->bc) {
489:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
490:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
491:   }
492:   return(0);
493: }

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

498:   Not collective

500:   Input Parameter:
501: . s - the PetscSection

503:   Output Parameters:
504: + pStart - the first point
505: - pEnd - one past the last point

507:   Level: intermediate

509: .seealso: PetscSectionSetChart(), PetscSectionCreate()
510: @*/
511: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
512: {
515:   if (pStart) *pStart = s->pStart;
516:   if (pEnd)   *pEnd   = s->pEnd;
517:   return(0);
518: }

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

523:   Not collective

525:   Input Parameters:
526: + s - the PetscSection
527: . pStart - the first point
528: - pEnd - one past the last point

530:   Level: intermediate

532: .seealso: PetscSectionGetChart(), PetscSectionCreate()
533: @*/
534: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
535: {
536:   PetscInt       f;

541:   /* Cannot Reset() because it destroys field information */
542:   s->setup = PETSC_FALSE;
543:   PetscSectionDestroy(&s->bc);
544:   PetscFree(s->bcIndices);
545:   PetscFree2(s->atlasDof, s->atlasOff);

547:   s->pStart = pStart;
548:   s->pEnd   = pEnd;
549:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
550:   PetscArrayzero(s->atlasDof, pEnd - pStart);
551:   for (f = 0; f < s->numFields; ++f) {
552:     PetscSectionSetChart(s->field[f], pStart, pEnd);
553:   }
554:   return(0);
555: }

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

560:   Not collective

562:   Input Parameter:
563: . s - the PetscSection

565:   Output Parameters:
566: . perm - The permutation as an IS

568:   Level: intermediate

570: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
571: @*/
572: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
573: {
577:   return(0);
578: }

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

583:   Not collective

585:   Input Parameters:
586: + s - the PetscSection
587: - perm - the permutation of points

589:   Level: intermediate

591: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
592: @*/
593: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
594: {

600:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
601:   if (s->perm != perm) {
602:     ISDestroy(&s->perm);
603:     if (perm) {
604:       s->perm = perm;
605:       PetscObjectReference((PetscObject) s->perm);
606:     }
607:   }
608:   return(0);
609: }

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

614:   Not collective

616:   Input Parameter:
617: . s - the PetscSection

619:   Output Parameter:
620: . pm - the flag for point major ordering

622:   Level: intermediate

624: .seealso: PetscSectionSetPointMajor()
625: @*/
626: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
627: {
631:   *pm = s->pointMajor;
632:   return(0);
633: }

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

638:   Not collective

640:   Input Parameters:
641: + s  - the PetscSection
642: - pm - the flag for point major ordering

644:   Not collective

646:   Level: intermediate

648: .seealso: PetscSectionGetPointMajor()
649: @*/
650: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
651: {
654:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
655:   s->pointMajor = pm;
656:   return(0);
657: }

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

662:   Not collective

664:   Input Parameters:
665: + s - the PetscSection
666: - point - the point

668:   Output Parameter:
669: . numDof - the number of dof

671:   Level: intermediate

673: .seealso: PetscSectionSetDof(), PetscSectionCreate()
674: @*/
675: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
676: {
680: #if defined(PETSC_USE_DEBUG)
681:   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);
682: #endif
683:   *numDof = s->atlasDof[point - s->pStart];
684:   return(0);
685: }

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

690:   Not collective

692:   Input Parameters:
693: + s - the PetscSection
694: . point - the point
695: - numDof - the number of dof

697:   Level: intermediate

699: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
700: @*/
701: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
702: {
705:   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);
706:   s->atlasDof[point - s->pStart] = numDof;
707:   return(0);
708: }

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

713:   Not collective

715:   Input Parameters:
716: + s - the PetscSection
717: . point - the point
718: - numDof - the number of additional dof

720:   Level: intermediate

722: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
723: @*/
724: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
725: {
728: #if defined(PETSC_USE_DEBUG)
729:   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);
730: #endif
731:   s->atlasDof[point - s->pStart] += numDof;
732:   return(0);
733: }

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

738:   Not collective

740:   Input Parameters:
741: + s - the PetscSection
742: . point - the point
743: - field - the field

745:   Output Parameter:
746: . numDof - the number of dof

748:   Level: intermediate

750: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
751: @*/
752: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
753: {

758:   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);
759:   PetscSectionGetDof(s->field[field], point, numDof);
760:   return(0);
761: }

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

766:   Not collective

768:   Input Parameters:
769: + s - the PetscSection
770: . point - the point
771: . field - the field
772: - numDof - the number of dof

774:   Level: intermediate

776: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
777: @*/
778: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
779: {

784:   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);
785:   PetscSectionSetDof(s->field[field], point, numDof);
786:   return(0);
787: }

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

792:   Not collective

794:   Input Parameters:
795: + s - the PetscSection
796: . point - the point
797: . field - the field
798: - numDof - the number of dof

800:   Level: intermediate

802: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
803: @*/
804: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
805: {

810:   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);
811:   PetscSectionAddDof(s->field[field], point, numDof);
812:   return(0);
813: }

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

818:   Not collective

820:   Input Parameters:
821: + s - the PetscSection
822: - point - the point

824:   Output Parameter:
825: . numDof - the number of dof which are fixed by constraints

827:   Level: intermediate

829: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
830: @*/
831: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
832: {

838:   if (s->bc) {
839:     PetscSectionGetDof(s->bc, point, numDof);
840:   } else *numDof = 0;
841:   return(0);
842: }

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

847:   Not collective

849:   Input Parameters:
850: + s - the PetscSection
851: . point - the point
852: - numDof - the number of dof which are fixed by constraints

854:   Level: intermediate

856: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
857: @*/
858: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
859: {

864:   if (numDof) {
865:     PetscSectionCheckConstraints_Static(s);
866:     PetscSectionSetDof(s->bc, point, numDof);
867:   }
868:   return(0);
869: }

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

874:   Not collective

876:   Input Parameters:
877: + s - the PetscSection
878: . point - the point
879: - numDof - the number of additional dof which are fixed by constraints

881:   Level: intermediate

883: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
884: @*/
885: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
886: {

891:   if (numDof) {
892:     PetscSectionCheckConstraints_Static(s);
893:     PetscSectionAddDof(s->bc, point, numDof);
894:   }
895:   return(0);
896: }

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

901:   Not collective

903:   Input Parameters:
904: + s - the PetscSection
905: . point - the point
906: - field - the field

908:   Output Parameter:
909: . numDof - the number of dof which are fixed by constraints

911:   Level: intermediate

913: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
914: @*/
915: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
916: {

922:   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);
923:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
924:   return(0);
925: }

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

930:   Not collective

932:   Input Parameters:
933: + s - the PetscSection
934: . point - the point
935: . field - the field
936: - numDof - the number of dof which are fixed by constraints

938:   Level: intermediate

940: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
941: @*/
942: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
943: {

948:   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);
949:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
950:   return(0);
951: }

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

956:   Not collective

958:   Input Parameters:
959: + s - the PetscSection
960: . point - the point
961: . field - the field
962: - numDof - the number of additional dof which are fixed by constraints

964:   Level: intermediate

966: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
967: @*/
968: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
969: {

974:   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);
975:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
976:   return(0);
977: }

979: /*@
980:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

982:   Not collective

984:   Input Parameter:
985: . s - the PetscSection

987:   Level: advanced

989: .seealso: PetscSectionSetUp(), PetscSectionCreate()
990: @*/
991: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
992: {

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

1000:     PetscSectionSetUp(s->bc);
1001:     PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1002:   }
1003:   return(0);
1004: }

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

1009:   Not collective

1011:   Input Parameter:
1012: . s - the PetscSection

1014:   Level: intermediate

1016: .seealso: PetscSectionCreate()
1017: @*/
1018: PetscErrorCode PetscSectionSetUp(PetscSection s)
1019: {
1020:   const PetscInt *pind   = NULL;
1021:   PetscInt        offset = 0, foff, p, f;
1022:   PetscErrorCode  ierr;

1026:   if (s->setup) return(0);
1027:   s->setup = PETSC_TRUE;
1028:   /* Set offsets and field offsets for all points */
1029:   /*   Assume that all fields have the same chart */
1030:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1031:   if (s->pointMajor) {
1032:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1033:       const PetscInt q = pind ? pind[p] : p;

1035:       /* Set point offset */
1036:       s->atlasOff[q] = offset;
1037:       offset        += s->atlasDof[q];
1038:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
1039:       /* Set field offset */
1040:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1041:         PetscSection sf = s->field[f];

1043:         sf->atlasOff[q] = foff;
1044:         foff += sf->atlasDof[q];
1045:       }
1046:     }
1047:   } else {
1048:     /* Set field offsets for all points */
1049:     for (f = 0; f < s->numFields; ++f) {
1050:       PetscSection sf = s->field[f];

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

1055:         sf->atlasOff[q] = offset;
1056:         offset += sf->atlasDof[q];
1057:       }
1058:     }
1059:     /* Disable point offsets since these are unused */
1060:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1061:       s->atlasOff[p] = -1;
1062:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
1063:     }
1064:   }
1065:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1066:   /* Setup BC sections */
1067:   PetscSectionSetUpBC(s);
1068:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1069:   return(0);
1070: }

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

1075:   Not collective

1077:   Input Parameters:
1078: . s - the PetscSection

1080:   Output Parameter:
1081: . maxDof - the maximum dof

1083:   Level: intermediate

1085: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1086: @*/
1087: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1088: {
1092:   *maxDof = s->maxDof;
1093:   return(0);
1094: }

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

1099:   Not collective

1101:   Input Parameter:
1102: . s - the PetscSection

1104:   Output Parameter:
1105: . size - the size of an array which can hold all the dofs

1107:   Level: intermediate

1109: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1110: @*/
1111: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1112: {
1113:   PetscInt p, n = 0;

1118:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1119:   *size = n;
1120:   return(0);
1121: }

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

1126:   Not collective

1128:   Input Parameters:
1129: + s - the PetscSection
1130: - point - the point

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

1135:   Level: intermediate

1137: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1138: @*/
1139: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1140: {
1141:   PetscInt p, n = 0;

1146:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1147:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1148:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1149:   }
1150:   *size = n;
1151:   return(0);
1152: }

1154: /*@
1155:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1156:   the local section and an SF describing the section point overlap.

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

1164:   Output Parameter:
1165: . gsection - The PetscSection for the global field layout

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

1169:   Level: intermediate

1171: .seealso: PetscSectionCreate()
1172: @*/
1173: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1174: {
1175:   PetscSection    gs;
1176:   const PetscInt *pind = NULL;
1177:   PetscInt       *recv = NULL, *neg = NULL;
1178:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1179:   PetscErrorCode  ierr;

1187:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1188:   PetscSectionGetChart(s, &pStart, &pEnd);
1189:   PetscSectionSetChart(gs, pStart, pEnd);
1190:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1191:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1192:   /* Must allocate for all points visible to SF, which may be more than this section */
1193:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1194:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1195:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1196:     if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1197:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1198:     PetscArrayzero(neg,nroots);
1199:   }
1200:   /* Mark all local points with negative dof */
1201:   for (p = pStart; p < pEnd; ++p) {
1202:     PetscSectionGetDof(s, p, &dof);
1203:     PetscSectionSetDof(gs, p, dof);
1204:     PetscSectionGetConstraintDof(s, p, &cdof);
1205:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1206:     if (neg) neg[p] = -(dof+1);
1207:   }
1208:   PetscSectionSetUpBC(gs);
1209:   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]);}
1210:   if (nroots >= 0) {
1211:     PetscArrayzero(recv,nlocal);
1212:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1213:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1214:     for (p = pStart; p < pEnd; ++p) {
1215:       if (recv[p] < 0) {
1216:         gs->atlasDof[p-pStart] = recv[p];
1217:         PetscSectionGetDof(s, p, &dof);
1218:         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);
1219:       }
1220:     }
1221:   }
1222:   /* Calculate new sizes, get process offset, and calculate point offsets */
1223:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1224:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1225:     const PetscInt q = pind ? pind[p] : p;

1227:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1228:     gs->atlasOff[q] = off;
1229:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1230:   }
1231:   if (!localOffsets) {
1232:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1233:     globalOff -= off;
1234:   }
1235:   for (p = pStart, off = 0; p < pEnd; ++p) {
1236:     gs->atlasOff[p-pStart] += globalOff;
1237:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1238:   }
1239:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1240:   /* Put in negative offsets for ghost points */
1241:   if (nroots >= 0) {
1242:     PetscArrayzero(recv,nlocal);
1243:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1244:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1245:     for (p = pStart; p < pEnd; ++p) {
1246:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1247:     }
1248:   }
1249:   PetscFree2(neg,recv);
1250:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1251:   *gsection = gs;
1252:   return(0);
1253: }

1255: /*@
1256:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1257:   the local section and an SF describing the section point overlap.

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

1266:   Output Parameter:
1267:   . gsection - The PetscSection for the global field layout

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

1271:   Level: advanced

1273: .seealso: PetscSectionCreate()
1274: @*/
1275: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1276: {
1277:   const PetscInt *pind = NULL;
1278:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1279:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1280:   PetscErrorCode  ierr;

1286:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1287:   PetscSectionGetChart(s, &pStart, &pEnd);
1288:   PetscSectionSetChart(*gsection, pStart, pEnd);
1289:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1290:   if (nroots >= 0) {
1291:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1292:     PetscCalloc1(nroots, &neg);
1293:     if (nroots > pEnd-pStart) {
1294:       PetscCalloc1(nroots, &tmpOff);
1295:     } else {
1296:       tmpOff = &(*gsection)->atlasDof[-pStart];
1297:     }
1298:   }
1299:   /* Mark ghost points with negative dof */
1300:   for (p = pStart; p < pEnd; ++p) {
1301:     for (e = 0; e < numExcludes; ++e) {
1302:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1303:         PetscSectionSetDof(*gsection, p, 0);
1304:         break;
1305:       }
1306:     }
1307:     if (e < numExcludes) continue;
1308:     PetscSectionGetDof(s, p, &dof);
1309:     PetscSectionSetDof(*gsection, p, dof);
1310:     PetscSectionGetConstraintDof(s, p, &cdof);
1311:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1312:     if (neg) neg[p] = -(dof+1);
1313:   }
1314:   PetscSectionSetUpBC(*gsection);
1315:   if (nroots >= 0) {
1316:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1317:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1318:     if (nroots > pEnd - pStart) {
1319:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1320:     }
1321:   }
1322:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1323:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1324:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1325:     const PetscInt q = pind ? pind[p] : p;

1327:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1328:     (*gsection)->atlasOff[q] = off;
1329:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1330:   }
1331:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1332:   globalOff -= off;
1333:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1334:     (*gsection)->atlasOff[p] += globalOff;
1335:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1336:   }
1337:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1338:   /* Put in negative offsets for ghost points */
1339:   if (nroots >= 0) {
1340:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1341:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1342:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1343:     if (nroots > pEnd - pStart) {
1344:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1345:     }
1346:   }
1347:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1348:   PetscFree(neg);
1349:   return(0);
1350: }

1352: /*@
1353:   PetscSectionGetPointLayout - Get the PetscLayout associated with the section points

1355:   Collective on comm

1357:   Input Parameters:
1358: + comm - The MPI_Comm
1359: - s    - The PetscSection

1361:   Output Parameter:
1362: . layout - The point layout for the section

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

1366:   Level: advanced

1368: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1369: @*/
1370: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1371: {
1372:   PetscInt       pStart, pEnd, p, localSize = 0;

1376:   PetscSectionGetChart(s, &pStart, &pEnd);
1377:   for (p = pStart; p < pEnd; ++p) {
1378:     PetscInt dof;

1380:     PetscSectionGetDof(s, p, &dof);
1381:     if (dof > 0) ++localSize;
1382:   }
1383:   PetscLayoutCreate(comm, layout);
1384:   PetscLayoutSetLocalSize(*layout, localSize);
1385:   PetscLayoutSetBlockSize(*layout, 1);
1386:   PetscLayoutSetUp(*layout);
1387:   return(0);
1388: }

1390: /*@
1391:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

1393:   Collective on comm

1395:   Input Parameters:
1396: + comm - The MPI_Comm
1397: - s    - The PetscSection

1399:   Output Parameter:
1400: . layout - The dof layout for the section

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

1404:   Level: advanced

1406: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1407: @*/
1408: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1409: {
1410:   PetscInt       pStart, pEnd, p, localSize = 0;

1416:   PetscSectionGetChart(s, &pStart, &pEnd);
1417:   for (p = pStart; p < pEnd; ++p) {
1418:     PetscInt dof,cdof;

1420:     PetscSectionGetDof(s, p, &dof);
1421:     PetscSectionGetConstraintDof(s, p, &cdof);
1422:     if (dof-cdof > 0) localSize += dof-cdof;
1423:   }
1424:   PetscLayoutCreate(comm, layout);
1425:   PetscLayoutSetLocalSize(*layout, localSize);
1426:   PetscLayoutSetBlockSize(*layout, 1);
1427:   PetscLayoutSetUp(*layout);
1428:   return(0);
1429: }

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

1434:   Not collective

1436:   Input Parameters:
1437: + s - the PetscSection
1438: - point - the point

1440:   Output Parameter:
1441: . offset - the offset

1443:   Level: intermediate

1445: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1446: @*/
1447: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1448: {
1452: #if defined(PETSC_USE_DEBUG)
1453:   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);
1454: #endif
1455:   *offset = s->atlasOff[point - s->pStart];
1456:   return(0);
1457: }

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

1462:   Not collective

1464:   Input Parameters:
1465: + s - the PetscSection
1466: . point - the point
1467: - offset - the offset

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

1471:   Level: intermediate

1473: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1474: @*/
1475: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1476: {
1479:   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);
1480:   s->atlasOff[point - s->pStart] = offset;
1481:   return(0);
1482: }

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

1487:   Not collective

1489:   Input Parameters:
1490: + s - the PetscSection
1491: . point - the point
1492: - field - the field

1494:   Output Parameter:
1495: . offset - the offset

1497:   Level: intermediate

1499: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1500: @*/
1501: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1502: {

1508:   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);
1509:   PetscSectionGetOffset(s->field[field], point, offset);
1510:   return(0);
1511: }

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

1516:   Not collective

1518:   Input Parameters:
1519: + s - the PetscSection
1520: . point - the point
1521: . field - the field
1522: - offset - the offset

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

1526:   Level: intermediate

1528: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1529: @*/
1530: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1531: {

1536:   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);
1537:   PetscSectionSetOffset(s->field[field], point, offset);
1538:   return(0);
1539: }

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

1544:   Not collective

1546:   Input Parameters:
1547: + s - the PetscSection
1548: . point - the point
1549: - field - the field

1551:   Output Parameter:
1552: . offset - the offset

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

1557:   Level: advanced

1559: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1560: @*/
1561: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1562: {
1563:   PetscInt       off, foff;

1569:   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);
1570:   PetscSectionGetOffset(s, point, &off);
1571:   PetscSectionGetOffset(s->field[field], point, &foff);
1572:   *offset = foff - off;
1573:   return(0);
1574: }

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

1579:   Not collective

1581:   Input Parameter:
1582: . s - the PetscSection

1584:   Output Parameters:
1585: + start - the minimum offset
1586: - end   - one more than the maximum offset

1588:   Level: intermediate

1590: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1591: @*/
1592: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1593: {
1594:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1599:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1600:   PetscSectionGetChart(s, &pStart, &pEnd);
1601:   for (p = 0; p < pEnd-pStart; ++p) {
1602:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1604:     if (off >= 0) {
1605:       os = PetscMin(os, off);
1606:       oe = PetscMax(oe, off+dof);
1607:     }
1608:   }
1609:   if (start) *start = os;
1610:   if (end)   *end   = oe;
1611:   return(0);
1612: }

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

1617:   Collective on s

1619:   Input Parameter:
1620: + s      - the PetscSection
1621: . len    - the number of subfields
1622: - fields - the subfield numbers

1624:   Output Parameter:
1625: . subs   - the subsection

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

1629:   Level: advanced

1631: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1632: @*/
1633: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1634: {
1635:   PetscInt       nF, f, pStart, pEnd, p, maxCdof = 0;

1639:   if (!len) return(0);
1643:   PetscSectionGetNumFields(s, &nF);
1644:   if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1645:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1646:   PetscSectionSetNumFields(*subs, len);
1647:   for (f = 0; f < len; ++f) {
1648:     const char *name   = NULL;
1649:     PetscInt   numComp = 0;

1651:     PetscSectionGetFieldName(s, fields[f], &name);
1652:     PetscSectionSetFieldName(*subs, f, name);
1653:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1654:     PetscSectionSetFieldComponents(*subs, f, numComp);
1655:   }
1656:   PetscSectionGetChart(s, &pStart, &pEnd);
1657:   PetscSectionSetChart(*subs, pStart, pEnd);
1658:   for (p = pStart; p < pEnd; ++p) {
1659:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1661:     for (f = 0; f < len; ++f) {
1662:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1663:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1664:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1665:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1666:       dof  += fdof;
1667:       cdof += cfdof;
1668:     }
1669:     PetscSectionSetDof(*subs, p, dof);
1670:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1671:     maxCdof = PetscMax(cdof, maxCdof);
1672:   }
1673:   PetscSectionSetUp(*subs);
1674:   if (maxCdof) {
1675:     PetscInt *indices;

1677:     PetscMalloc1(maxCdof, &indices);
1678:     for (p = pStart; p < pEnd; ++p) {
1679:       PetscInt cdof;

1681:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1682:       if (cdof) {
1683:         const PetscInt *oldIndices = NULL;
1684:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1689:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1690:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1691:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1692:           /* This can be sped up if we assume sorted fields */
1693:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1694:             PetscInt oldfdof = 0;
1695:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1696:             oldFoff += oldfdof;
1697:           }
1698:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1699:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1700:           numConst += cfdof;
1701:           fOff     += fdof;
1702:         }
1703:         PetscSectionSetConstraintIndices(*subs, p, indices);
1704:       }
1705:     }
1706:     PetscFree(indices);
1707:   }
1708:   return(0);
1709: }

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

1714:   Collective on s

1716:   Input Parameters:
1717: + s     - the input sections
1718: - len   - the number of input sections

1720:   Output Parameter:
1721: . supers - the supersection

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

1725:   Level: advanced

1727: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1728: @*/
1729: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1730: {
1731:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1735:   if (!len) return(0);
1736:   for (i = 0; i < len; ++i) {
1737:     PetscInt nf, pStarti, pEndi;

1739:     PetscSectionGetNumFields(s[i], &nf);
1740:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1741:     pStart = PetscMin(pStart, pStarti);
1742:     pEnd   = PetscMax(pEnd,   pEndi);
1743:     Nf += nf;
1744:   }
1745:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1746:   PetscSectionSetNumFields(*supers, Nf);
1747:   for (i = 0, f = 0; i < len; ++i) {
1748:     PetscInt nf, fi;

1750:     PetscSectionGetNumFields(s[i], &nf);
1751:     for (fi = 0; fi < nf; ++fi, ++f) {
1752:       const char *name   = NULL;
1753:       PetscInt   numComp = 0;

1755:       PetscSectionGetFieldName(s[i], fi, &name);
1756:       PetscSectionSetFieldName(*supers, f, name);
1757:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1758:       PetscSectionSetFieldComponents(*supers, f, numComp);
1759:     }
1760:   }
1761:   PetscSectionSetChart(*supers, pStart, pEnd);
1762:   for (p = pStart; p < pEnd; ++p) {
1763:     PetscInt dof = 0, cdof = 0;

1765:     for (i = 0, f = 0; i < len; ++i) {
1766:       PetscInt nf, fi, pStarti, pEndi;
1767:       PetscInt fdof = 0, cfdof = 0;

1769:       PetscSectionGetNumFields(s[i], &nf);
1770:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1771:       if ((p < pStarti) || (p >= pEndi)) continue;
1772:       for (fi = 0; fi < nf; ++fi, ++f) {
1773:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1774:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1775:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1776:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1777:         dof  += fdof;
1778:         cdof += cfdof;
1779:       }
1780:     }
1781:     PetscSectionSetDof(*supers, p, dof);
1782:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1783:     maxCdof = PetscMax(cdof, maxCdof);
1784:   }
1785:   PetscSectionSetUp(*supers);
1786:   if (maxCdof) {
1787:     PetscInt *indices;

1789:     PetscMalloc1(maxCdof, &indices);
1790:     for (p = pStart; p < pEnd; ++p) {
1791:       PetscInt cdof;

1793:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1794:       if (cdof) {
1795:         PetscInt dof, numConst = 0, fOff = 0;

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

1801:           PetscSectionGetNumFields(s[i], &nf);
1802:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1803:           if ((p < pStarti) || (p >= pEndi)) continue;
1804:           for (fi = 0; fi < nf; ++fi, ++f) {
1805:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1806:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1807:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1808:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1809:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1810:             numConst += cfdof;
1811:           }
1812:           PetscSectionGetDof(s[i], p, &dof);
1813:           fOff += dof;
1814:         }
1815:         PetscSectionSetConstraintIndices(*supers, p, indices);
1816:       }
1817:     }
1818:     PetscFree(indices);
1819:   }
1820:   return(0);
1821: }

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

1826:   Collective on s

1828:   Input Parameters:
1829: + s           - the PetscSection
1830: - subpointMap - a sorted list of points in the original mesh which are in the submesh

1832:   Output Parameter:
1833: . subs - the subsection

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

1837:   Level: advanced

1839: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1840: @*/
1841: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1842: {
1843:   const PetscInt *points = NULL, *indices = NULL;
1844:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1851:   PetscSectionGetNumFields(s, &numFields);
1852:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1853:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1854:   for (f = 0; f < numFields; ++f) {
1855:     const char *name   = NULL;
1856:     PetscInt   numComp = 0;

1858:     PetscSectionGetFieldName(s, f, &name);
1859:     PetscSectionSetFieldName(*subs, f, name);
1860:     PetscSectionGetFieldComponents(s, f, &numComp);
1861:     PetscSectionSetFieldComponents(*subs, f, numComp);
1862:   }
1863:   /* For right now, we do not try to squeeze the subchart */
1864:   if (subpointMap) {
1865:     ISGetSize(subpointMap, &numSubpoints);
1866:     ISGetIndices(subpointMap, &points);
1867:   }
1868:   PetscSectionGetChart(s, &pStart, &pEnd);
1869:   PetscSectionSetChart(*subs, 0, numSubpoints);
1870:   for (p = pStart; p < pEnd; ++p) {
1871:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1873:     PetscFindInt(p, numSubpoints, points, &subp);
1874:     if (subp < 0) continue;
1875:     for (f = 0; f < numFields; ++f) {
1876:       PetscSectionGetFieldDof(s, p, f, &fdof);
1877:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1878:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1879:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1880:     }
1881:     PetscSectionGetDof(s, p, &dof);
1882:     PetscSectionSetDof(*subs, subp, dof);
1883:     PetscSectionGetConstraintDof(s, p, &cdof);
1884:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1885:   }
1886:   PetscSectionSetUp(*subs);
1887:   /* Change offsets to original offsets */
1888:   for (p = pStart; p < pEnd; ++p) {
1889:     PetscInt off, foff = 0;

1891:     PetscFindInt(p, numSubpoints, points, &subp);
1892:     if (subp < 0) continue;
1893:     for (f = 0; f < numFields; ++f) {
1894:       PetscSectionGetFieldOffset(s, p, f, &foff);
1895:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1896:     }
1897:     PetscSectionGetOffset(s, p, &off);
1898:     PetscSectionSetOffset(*subs, subp, off);
1899:   }
1900:   /* Copy constraint indices */
1901:   for (subp = 0; subp < numSubpoints; ++subp) {
1902:     PetscInt cdof;

1904:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1905:     if (cdof) {
1906:       for (f = 0; f < numFields; ++f) {
1907:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1908:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1909:       }
1910:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1911:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1912:     }
1913:   }
1914:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1915:   return(0);
1916: }

1918: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1919: {
1920:   PetscInt       p;
1921:   PetscMPIInt    rank;

1925:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1926:   PetscViewerASCIIPushSynchronized(viewer);
1927:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1928:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1929:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1930:       PetscInt b;

1932:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1933:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1934:         PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
1935:       }
1936:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1937:     } else {
1938:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1939:     }
1940:   }
1941:   PetscViewerFlush(viewer);
1942:   PetscViewerASCIIPopSynchronized(viewer);
1943:   if (s->sym) {
1944:     PetscViewerASCIIPushTab(viewer);
1945:     PetscSectionSymView(s->sym,viewer);
1946:     PetscViewerASCIIPopTab(viewer);
1947:   }
1948:   return(0);
1949: }

1951: /*@C
1952:    PetscSectionViewFromOptions - View from Options

1954:    Collective on PetscSection

1956:    Input Parameters:
1957: +  A - the PetscSection object to view
1958: .  obj - Optional object
1959: -  name - command line option

1961:    Level: intermediate
1962: .seealso:  PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
1963: @*/
1964: PetscErrorCode  PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
1965: {

1970:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
1971:   return(0);
1972: }

1974: /*@C
1975:   PetscSectionView - Views a PetscSection

1977:   Collective on PetscSection

1979:   Input Parameters:
1980: + s - the PetscSection object to view
1981: - v - the viewer

1983:   Level: beginner

1985: .seealso PetscSectionCreate(), PetscSectionDestroy()
1986: @*/
1987: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1988: {
1989:   PetscBool      isascii;
1990:   PetscInt       f;

1995:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1997:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1998:   if (isascii) {
1999:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2000:     if (s->numFields) {
2001:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2002:       for (f = 0; f < s->numFields; ++f) {
2003:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
2004:         PetscSectionView_ASCII(s->field[f], viewer);
2005:       }
2006:     } else {
2007:       PetscSectionView_ASCII(s, viewer);
2008:     }
2009:   }
2010:   return(0);
2011: }

2013: /*@
2014:   PetscSectionReset - Frees all section data.

2016:   Not collective

2018:   Input Parameters:
2019: . s - the PetscSection

2021:   Level: beginner

2023: .seealso: PetscSection, PetscSectionCreate()
2024: @*/
2025: PetscErrorCode PetscSectionReset(PetscSection s)
2026: {
2027:   PetscInt       f;

2032:   PetscFree(s->numFieldComponents);
2033:   for (f = 0; f < s->numFields; ++f) {
2034:     PetscSectionDestroy(&s->field[f]);
2035:     PetscFree(s->fieldNames[f]);
2036:   }
2037:   PetscFree(s->fieldNames);
2038:   PetscFree(s->field);
2039:   PetscSectionDestroy(&s->bc);
2040:   PetscFree(s->bcIndices);
2041:   PetscFree2(s->atlasDof, s->atlasOff);
2042:   PetscSectionDestroy(&s->clSection);
2043:   ISDestroy(&s->clPoints);
2044:   ISDestroy(&s->perm);
2045:   PetscFree(s->clPerm);
2046:   PetscFree(s->clInvPerm);
2047:   PetscSectionSymDestroy(&s->sym);
2048:   PetscSectionDestroy(&s->clSection);
2049:   ISDestroy(&s->clPoints);

2051:   s->pStart    = -1;
2052:   s->pEnd      = -1;
2053:   s->maxDof    = 0;
2054:   s->setup     = PETSC_FALSE;
2055:   s->numFields = 0;
2056:   s->clObj     = NULL;
2057:   return(0);
2058: }

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

2063:   Not collective

2065:   Input Parameters:
2066: . s - the PetscSection

2068:   Level: beginner

2070: .seealso: PetscSection, PetscSectionCreate()
2071: @*/
2072: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2073: {

2077:   if (!*s) return(0);
2079:   if (--((PetscObject)(*s))->refct > 0) {
2080:     *s = NULL;
2081:     return(0);
2082:   }
2083:   PetscSectionReset(*s);
2084:   PetscHeaderDestroy(s);
2085:   return(0);
2086: }

2088: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2089: {
2090:   const PetscInt p = point - s->pStart;

2094:   *values = &baseArray[s->atlasOff[p]];
2095:   return(0);
2096: }

2098: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2099: {
2100:   PetscInt       *array;
2101:   const PetscInt p           = point - s->pStart;
2102:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2103:   PetscInt       cDim        = 0;

2108:   PetscSectionGetConstraintDof(s, p, &cDim);
2109:   array = &baseArray[s->atlasOff[p]];
2110:   if (!cDim) {
2111:     if (orientation >= 0) {
2112:       const PetscInt dim = s->atlasDof[p];
2113:       PetscInt       i;

2115:       if (mode == INSERT_VALUES) {
2116:         for (i = 0; i < dim; ++i) array[i] = values[i];
2117:       } else {
2118:         for (i = 0; i < dim; ++i) array[i] += values[i];
2119:       }
2120:     } else {
2121:       PetscInt offset = 0;
2122:       PetscInt j      = -1, field, i;

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

2127:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2128:         offset += dim;
2129:       }
2130:     }
2131:   } else {
2132:     if (orientation >= 0) {
2133:       const PetscInt dim  = s->atlasDof[p];
2134:       PetscInt       cInd = 0, i;
2135:       const PetscInt *cDof;

2137:       PetscSectionGetConstraintIndices(s, point, &cDof);
2138:       if (mode == INSERT_VALUES) {
2139:         for (i = 0; i < dim; ++i) {
2140:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2141:           array[i] = values[i];
2142:         }
2143:       } else {
2144:         for (i = 0; i < dim; ++i) {
2145:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2146:           array[i] += values[i];
2147:         }
2148:       }
2149:     } else {
2150:       const PetscInt *cDof;
2151:       PetscInt       offset  = 0;
2152:       PetscInt       cOffset = 0;
2153:       PetscInt       j       = 0, field;

2155:       PetscSectionGetConstraintIndices(s, point, &cDof);
2156:       for (field = 0; field < s->numFields; ++field) {
2157:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2158:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2159:         const PetscInt sDim = dim - tDim;
2160:         PetscInt       cInd = 0, i,k;

2162:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2163:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2164:           array[j] = values[k];
2165:         }
2166:         offset  += dim;
2167:         cOffset += dim - tDim;
2168:       }
2169:     }
2170:   }
2171:   return(0);
2172: }

2174: /*@C
2175:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2177:   Not collective

2179:   Input Parameter:
2180: . s - The PetscSection

2182:   Output Parameter:
2183: . hasConstraints - flag indicating that the section has constrained dofs

2185:   Level: intermediate

2187: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2188: @*/
2189: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2190: {
2194:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2195:   return(0);
2196: }

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

2201:   Not collective

2203:   Input Parameters:
2204: + s     - The PetscSection
2205: - point - The point

2207:   Output Parameter:
2208: . indices - The constrained dofs

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

2212:   Level: intermediate

2214: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2215: @*/
2216: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2217: {

2222:   if (s->bc) {
2223:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2224:   } else *indices = NULL;
2225:   return(0);
2226: }

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

2231:   Not collective

2233:   Input Parameters:
2234: + s     - The PetscSection
2235: . point - The point
2236: - indices - The constrained dofs

2238:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2240:   Level: intermediate

2242: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2243: @*/
2244: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2245: {

2250:   if (s->bc) {
2251:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2252:   }
2253:   return(0);
2254: }

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

2259:   Not collective

2261:   Input Parameters:
2262: + s     - The PetscSection
2263: . field  - The field number
2264: - point - The point

2266:   Output Parameter:
2267: . indices - The constrained dofs

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

2271:   Level: intermediate

2273: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2274: @*/
2275: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2276: {

2281:   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);
2282:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2283:   return(0);
2284: }

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

2289:   Not collective

2291:   Input Parameters:
2292: + s       - The PetscSection
2293: . point   - The point
2294: . field   - The field number
2295: - indices - The constrained dofs

2297:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2299:   Level: intermediate

2301: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2302: @*/
2303: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2304: {

2309:   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);
2310:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2311:   return(0);
2312: }

2314: /*@
2315:   PetscSectionPermute - Reorder the section according to the input point permutation

2317:   Collective on PetscSection

2319:   Input Parameter:
2320: + section - The PetscSection object
2321: - perm - The point permutation, old point p becomes new point perm[p]

2323:   Output Parameter:
2324: . sectionNew - The permuted PetscSection

2326:   Level: intermediate

2328: .seealso: MatPermute()
2329: @*/
2330: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2331: {
2332:   PetscSection    s = section, sNew;
2333:   const PetscInt *perm;
2334:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
2335:   PetscErrorCode  ierr;

2341:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2342:   PetscSectionGetNumFields(s, &numFields);
2343:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2344:   for (f = 0; f < numFields; ++f) {
2345:     const char *name;
2346:     PetscInt    numComp;

2348:     PetscSectionGetFieldName(s, f, &name);
2349:     PetscSectionSetFieldName(sNew, f, name);
2350:     PetscSectionGetFieldComponents(s, f, &numComp);
2351:     PetscSectionSetFieldComponents(sNew, f, numComp);
2352:   }
2353:   ISGetLocalSize(permutation, &numPoints);
2354:   ISGetIndices(permutation, &perm);
2355:   PetscSectionGetChart(s, &pStart, &pEnd);
2356:   PetscSectionSetChart(sNew, pStart, pEnd);
2357:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2358:   for (p = pStart; p < pEnd; ++p) {
2359:     PetscInt dof, cdof;

2361:     PetscSectionGetDof(s, p, &dof);
2362:     PetscSectionSetDof(sNew, perm[p], dof);
2363:     PetscSectionGetConstraintDof(s, p, &cdof);
2364:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2365:     for (f = 0; f < numFields; ++f) {
2366:       PetscSectionGetFieldDof(s, p, f, &dof);
2367:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2368:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2369:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2370:     }
2371:   }
2372:   PetscSectionSetUp(sNew);
2373:   for (p = pStart; p < pEnd; ++p) {
2374:     const PetscInt *cind;
2375:     PetscInt        cdof;

2377:     PetscSectionGetConstraintDof(s, p, &cdof);
2378:     if (cdof) {
2379:       PetscSectionGetConstraintIndices(s, p, &cind);
2380:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2381:     }
2382:     for (f = 0; f < numFields; ++f) {
2383:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2384:       if (cdof) {
2385:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2386:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2387:       }
2388:     }
2389:   }
2390:   ISRestoreIndices(permutation, &perm);
2391:   *sectionNew = sNew;
2392:   return(0);
2393: }

2395: /* TODO: the next three functions should be moved to sf/utils */
2396:  #include <petsc/private/sfimpl.h>

2398: /*@C
2399:   PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF

2401:   Collective on sf

2403:   Input Parameters:
2404: + sf - The SF
2405: - rootSection - Section defined on root space

2407:   Output Parameters:
2408: + remoteOffsets - root offsets in leaf storage, or NULL
2409: - leafSection - Section defined on the leaf space

2411:   Level: advanced

2413: .seealso: PetscSFCreate()
2414: @*/
2415: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2416: {
2417:   PetscSF        embedSF;
2418:   const PetscInt *indices;
2419:   IS             selected;
2420:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f;
2421:   PetscBool      *sub, hasc;

2425:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2426:   PetscSectionGetNumFields(rootSection, &numFields);
2427:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2428:   PetscMalloc1(numFields+2, &sub);
2429:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2430:   for (f = 0; f < numFields; ++f) {
2431:     PetscSectionSym sym;
2432:     const char      *name   = NULL;
2433:     PetscInt        numComp = 0;

2435:     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2436:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2437:     PetscSectionGetFieldName(rootSection, f, &name);
2438:     PetscSectionGetFieldSym(rootSection, f, &sym);
2439:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2440:     PetscSectionSetFieldName(leafSection, f, name);
2441:     PetscSectionSetFieldSym(leafSection, f, sym);
2442:   }
2443:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2444:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2445:   rpEnd = PetscMin(rpEnd,nroots);
2446:   rpEnd = PetscMax(rpStart,rpEnd);
2447:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2448:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2449:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2450:   if (sub[0]) {
2451:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2452:     ISGetIndices(selected, &indices);
2453:     PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2454:     ISRestoreIndices(selected, &indices);
2455:     ISDestroy(&selected);
2456:   } else {
2457:     PetscObjectReference((PetscObject)sf);
2458:     embedSF = sf;
2459:   }
2460:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2461:   lpEnd++;

2463:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

2465:   /* Constrained dof section */
2466:   hasc = sub[1];
2467:   for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);

2469:   /* Could fuse these at the cost of copies and extra allocation */
2470:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2471:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2472:   if (sub[1]) {
2473:     PetscSectionCheckConstraints_Static(rootSection);
2474:     PetscSectionCheckConstraints_Static(leafSection);
2475:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2476:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2477:   }
2478:   for (f = 0; f < numFields; ++f) {
2479:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2480:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2481:     if (sub[2+f]) {
2482:       PetscSectionCheckConstraints_Static(rootSection->field[f]);
2483:       PetscSectionCheckConstraints_Static(leafSection->field[f]);
2484:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2485:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2486:     }
2487:   }
2488:   if (remoteOffsets) {
2489:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2490:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2491:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2492:   }
2493:   PetscSectionSetUp(leafSection);
2494:   if (hasc) { /* need to communicate bcIndices */
2495:     PetscSF  bcSF;
2496:     PetscInt *rOffBc;

2498:     PetscMalloc1(lpEnd - lpStart, &rOffBc);
2499:     if (sub[1]) {
2500:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2501:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2502:       PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2503:       PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2504:       PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2505:       PetscSFDestroy(&bcSF);
2506:     }
2507:     for (f = 0; f < numFields; ++f) {
2508:       if (sub[2+f]) {
2509:         PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2510:         PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2511:         PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2512:         PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2513:         PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2514:         PetscSFDestroy(&bcSF);
2515:       }
2516:     }
2517:     PetscFree(rOffBc);
2518:   }
2519:   PetscSFDestroy(&embedSF);
2520:   PetscFree(sub);
2521:   PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2522:   return(0);
2523: }

2525: /*@C
2526:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

2528:   Collective on sf

2530:   Input Parameters:
2531: + sf          - The SF
2532: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2533: - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)

2535:   Output Parameter:
2536: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL

2538:   Level: developer

2540: .seealso: PetscSFCreate()
2541: @*/
2542: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2543: {
2544:   PetscSF         embedSF;
2545:   const PetscInt *indices;
2546:   IS              selected;
2547:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2548:   PetscErrorCode  ierr;

2551:   *remoteOffsets = NULL;
2552:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2553:   if (numRoots < 0) return(0);
2554:   PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2555:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2556:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2557:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2558:   ISGetIndices(selected, &indices);
2559:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2560:   ISRestoreIndices(selected, &indices);
2561:   ISDestroy(&selected);
2562:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2563:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2564:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2565:   PetscSFDestroy(&embedSF);
2566:   PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2567:   return(0);
2568: }

2570: /*@C
2571:   PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points

2573:   Collective on sf

2575:   Input Parameters:
2576: + sf - The SF
2577: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2578: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2579: - leafSection - Data layout of local points for incoming data  (this is the distributed section)

2581:   Output Parameters:
2582: - sectionSF - The new SF

2584:   Note: Either rootSection or remoteOffsets can be specified

2586:   Level: advanced

2588: .seealso: PetscSFCreate()
2589: @*/
2590: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2591: {
2592:   MPI_Comm          comm;
2593:   const PetscInt    *localPoints;
2594:   const PetscSFNode *remotePoints;
2595:   PetscInt          lpStart, lpEnd;
2596:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2597:   PetscInt          *localIndices;
2598:   PetscSFNode       *remoteIndices;
2599:   PetscInt          i, ind;
2600:   PetscErrorCode    ierr;

2608:   PetscObjectGetComm((PetscObject)sf,&comm);
2609:   PetscSFCreate(comm, sectionSF);
2610:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2611:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2612:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2613:   if (numRoots < 0) return(0);
2614:   PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2615:   for (i = 0; i < numPoints; ++i) {
2616:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2617:     PetscInt dof;

2619:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2620:       PetscSectionGetDof(leafSection, localPoint, &dof);
2621:       numIndices += dof;
2622:     }
2623:   }
2624:   PetscMalloc1(numIndices, &localIndices);
2625:   PetscMalloc1(numIndices, &remoteIndices);
2626:   /* Create new index graph */
2627:   for (i = 0, ind = 0; i < numPoints; ++i) {
2628:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2629:     PetscInt rank       = remotePoints[i].rank;

2631:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2632:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2633:       PetscInt loff, dof, d;

2635:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2636:       PetscSectionGetDof(leafSection, localPoint, &dof);
2637:       for (d = 0; d < dof; ++d, ++ind) {
2638:         localIndices[ind]        = loff+d;
2639:         remoteIndices[ind].rank  = rank;
2640:         remoteIndices[ind].index = remoteOffset+d;
2641:       }
2642:     }
2643:   }
2644:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2645:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2646:   PetscSFSetUp(*sectionSF);
2647:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2648:   return(0);
2649: }

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

2654:   Collective on section

2656:   Input Parameters:
2657: + section   - The PetscSection
2658: . obj       - A PetscObject which serves as the key for this index
2659: . clSection - Section giving the size of the closure of each point
2660: - clPoints  - IS giving the points in each closure

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

2664:   Level: advanced

2666: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2667: @*/
2668: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2669: {

2676:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2677:   section->clObj     = obj;
2678:   PetscObjectReference((PetscObject)clSection);
2679:   PetscObjectReference((PetscObject)clPoints);
2680:   PetscSectionDestroy(&section->clSection);
2681:   ISDestroy(&section->clPoints);
2682:   section->clSection = clSection;
2683:   section->clPoints  = clPoints;
2684:   return(0);
2685: }

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

2690:   Collective on section

2692:   Input Parameters:
2693: + section   - The PetscSection
2694: - obj       - A PetscObject which serves as the key for this index

2696:   Output Parameters:
2697: + clSection - Section giving the size of the closure of each point
2698: - clPoints  - IS giving the points in each closure

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

2702:   Level: advanced

2704: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2705: @*/
2706: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2707: {
2709:   if (section->clObj == obj) {
2710:     if (clSection) *clSection = section->clSection;
2711:     if (clPoints)  *clPoints  = section->clPoints;
2712:   } else {
2713:     if (clSection) *clSection = NULL;
2714:     if (clPoints)  *clPoints  = NULL;
2715:   }
2716:   return(0);
2717: }

2719: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2720: {
2721:   PetscInt       i;

2725:   if (section->clObj != obj) {
2726:     PetscSectionDestroy(&section->clSection);
2727:     ISDestroy(&section->clPoints);
2728:   }
2729:   section->clObj  = obj;
2730:   PetscFree(section->clPerm);
2731:   PetscFree(section->clInvPerm);
2732:   section->clSize = clSize;
2733:   if (mode == PETSC_COPY_VALUES) {
2734:     PetscMalloc1(clSize, &section->clPerm);
2735:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2736:     PetscArraycpy(section->clPerm, clPerm, clSize);
2737:   } else if (mode == PETSC_OWN_POINTER) {
2738:     section->clPerm = clPerm;
2739:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2740:   PetscMalloc1(clSize, &section->clInvPerm);
2741:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2742:   return(0);
2743: }

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

2748:   Not Collective

2750:   Input Parameters:
2751: + section - The PetscSection
2752: . obj     - A PetscObject which serves as the key for this index
2753: - perm    - Permutation of the cell dof closure

2755:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2756:   other points (like faces).

2758:   Level: intermediate

2760: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2761: @*/
2762: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2763: {
2764:   const PetscInt *clPerm = NULL;
2765:   PetscInt        clSize = 0;
2766:   PetscErrorCode  ierr;

2769:   if (perm) {
2770:     ISGetLocalSize(perm, &clSize);
2771:     ISGetIndices(perm, &clPerm);
2772:   }
2773:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2774:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2775:   return(0);
2776: }

2778: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2779: {
2781:   if (section->clObj == obj) {
2782:     if (size) *size = section->clSize;
2783:     if (perm) *perm = section->clPerm;
2784:   } else {
2785:     if (size) *size = 0;
2786:     if (perm) *perm = NULL;
2787:   }
2788:   return(0);
2789: }

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

2794:   Not collective

2796:   Input Parameters:
2797: + section   - The PetscSection
2798: - obj       - A PetscObject which serves as the key for this index

2800:   Output Parameter:
2801: . perm - The dof closure permutation

2803:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2804:   other points (like faces).

2806:   The user must destroy the IS that is returned.

2808:   Level: intermediate

2810: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2811: @*/
2812: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2813: {
2814:   const PetscInt *clPerm;
2815:   PetscInt        clSize;
2816:   PetscErrorCode  ierr;

2819:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2820:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2821:   return(0);
2822: }

2824: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2825: {
2827:   if (section->clObj == obj) {
2828:     if (size) *size = section->clSize;
2829:     if (perm) *perm = section->clInvPerm;
2830:   } else {
2831:     if (size) *size = 0;
2832:     if (perm) *perm = NULL;
2833:   }
2834:   return(0);
2835: }

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

2840:   Not collective

2842:   Input Parameters:
2843: + section   - The PetscSection
2844: - obj       - A PetscObject which serves as the key for this index

2846:   Output Parameters:
2847: + size - The dof closure size
2848: - perm - The dof closure permutation

2850:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2851:   other points (like faces).

2853:   The user must destroy the IS that is returned.

2855:   Level: intermediate

2857: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2858: @*/
2859: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2860: {
2861:   const PetscInt *clPerm;
2862:   PetscInt        clSize;
2863:   PetscErrorCode  ierr;

2866:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2867:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2868:   return(0);
2869: }

2871: /*@
2872:   PetscSectionGetField - Get the subsection associated with a single field

2874:   Input Parameters:
2875: + s     - The PetscSection
2876: - field - The field number

2878:   Output Parameter:
2879: . subs  - The subsection for the given field

2881:   Level: intermediate

2883: .seealso: PetscSectionSetNumFields()
2884: @*/
2885: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2886: {
2890:   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);
2891:   *subs = s->field[field];
2892:   return(0);
2893: }

2895: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2896: PetscFunctionList PetscSectionSymList = NULL;

2898: /*@
2899:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2901:   Collective

2903:   Input Parameter:
2904: . comm - the MPI communicator

2906:   Output Parameter:
2907: . sym - pointer to the new set of symmetries

2909:   Level: developer

2911: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2912: @*/
2913: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2914: {

2919:   ISInitializePackage();
2920:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2921:   return(0);
2922: }

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

2927:   Collective on PetscSectionSym

2929:   Input Parameters:
2930: + sym    - The section symmetry object
2931: - method - The name of the section symmetry type

2933:   Level: developer

2935: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2936: @*/
2937: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2938: {
2939:   PetscErrorCode (*r)(PetscSectionSym);
2940:   PetscBool      match;

2945:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2946:   if (match) return(0);

2948:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2949:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2950:   if (sym->ops->destroy) {
2951:     (*sym->ops->destroy)(sym);
2952:     sym->ops->destroy = NULL;
2953:   }
2954:   (*r)(sym);
2955:   PetscObjectChangeTypeName((PetscObject)sym,method);
2956:   return(0);
2957: }


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

2963:   Not Collective

2965:   Input Parameter:
2966: . sym  - The section symmetry

2968:   Output Parameter:
2969: . type - The index set type name

2971:   Level: developer

2973: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2974: @*/
2975: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2976: {
2980:   *type = ((PetscObject)sym)->type_name;
2981:   return(0);
2982: }

2984: /*@C
2985:   PetscSectionSymRegister - Adds a new section symmetry implementation

2987:   Not Collective

2989:   Input Parameters:
2990: + name        - The name of a new user-defined creation routine
2991: - create_func - The creation routine itself

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

2996:   Level: developer

2998: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2999: @*/
3000: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3001: {

3005:   ISInitializePackage();
3006:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3007:   return(0);
3008: }

3010: /*@
3011:    PetscSectionSymDestroy - Destroys a section symmetry.

3013:    Collective on PetscSectionSym

3015:    Input Parameters:
3016: .  sym - the section symmetry

3018:    Level: developer

3020: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3021: @*/
3022: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3023: {
3024:   SymWorkLink    link,next;

3028:   if (!*sym) return(0);
3030:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
3031:   if ((*sym)->ops->destroy) {
3032:     (*(*sym)->ops->destroy)(*sym);
3033:   }
3034:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3035:   for (link=(*sym)->workin; link; link=next) {
3036:     next = link->next;
3037:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3038:     PetscFree(link);
3039:   }
3040:   (*sym)->workin = NULL;
3041:   PetscHeaderDestroy(sym);
3042:   return(0);
3043: }

3045: /*@C
3046:    PetscSectionSymView - Displays a section symmetry

3048:    Collective on PetscSectionSym

3050:    Input Parameters:
3051: +  sym - the index set
3052: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

3054:    Level: developer

3056: .seealso: PetscViewerASCIIOpen()
3057: @*/
3058: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3059: {

3064:   if (!viewer) {
3065:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3066:   }
3069:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3070:   if (sym->ops->view) {
3071:     (*sym->ops->view)(sym,viewer);
3072:   }
3073:   return(0);
3074: }

3076: /*@
3077:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3079:   Collective on PetscSection

3081:   Input Parameters:
3082: + section - the section describing data layout
3083: - sym - the symmetry describing the affect of orientation on the access of the data

3085:   Level: developer

3087: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3088: @*/
3089: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3090: {

3095:   PetscSectionSymDestroy(&(section->sym));
3096:   if (sym) {
3099:     PetscObjectReference((PetscObject) sym);
3100:   }
3101:   section->sym = sym;
3102:   return(0);
3103: }

3105: /*@
3106:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3108:   Not collective

3110:   Input Parameters:
3111: . section - the section describing data layout

3113:   Output Parameters:
3114: . sym - the symmetry describing the affect of orientation on the access of the data

3116:   Level: developer

3118: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3119: @*/
3120: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3121: {
3124:   *sym = section->sym;
3125:   return(0);
3126: }

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

3131:   Collective on PetscSection

3133:   Input Parameters:
3134: + section - the section describing data layout
3135: . field - the field number
3136: - sym - the symmetry describing the affect of orientation on the access of the data

3138:   Level: developer

3140: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3141: @*/
3142: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3143: {

3148:   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);
3149:   PetscSectionSetSym(section->field[field],sym);
3150:   return(0);
3151: }

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

3156:   Collective on PetscSection

3158:   Input Parameters:
3159: + section - the section describing data layout
3160: - field - the field number

3162:   Output Parameters:
3163: . sym - the symmetry describing the affect of orientation on the access of the data

3165:   Level: developer

3167: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3168: @*/
3169: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3170: {
3173:   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);
3174:   *sym = section->field[field]->sym;
3175:   return(0);
3176: }

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

3181:   Not collective

3183:   Input Parameters:
3184: + section - the section
3185: . numPoints - the number of points
3186: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3187:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3188:     context, see DMPlexGetConeOrientation()).

3190:   Output Parameter:
3191: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3192: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3193:     identity).

3195:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3196: .vb
3197:      const PetscInt    **perms;
3198:      const PetscScalar **rots;
3199:      PetscInt            lOffset;

3201:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3202:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3203:        PetscInt           point = points[2*i], dof, sOffset;
3204:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3205:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3207:        PetscSectionGetDof(section,point,&dof);
3208:        PetscSectionGetOffset(section,point,&sOffset);

3210:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3211:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3212:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3213:        lOffset += dof;
3214:      }
3215:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3216: .ve

3218:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3219: .vb
3220:      const PetscInt    **perms;
3221:      const PetscScalar **rots;
3222:      PetscInt            lOffset;

3224:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3225:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3226:        PetscInt           point = points[2*i], dof, sOffset;
3227:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3228:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3230:        PetscSectionGetDof(section,point,&dof);
3231:        PetscSectionGetOffset(section,point,&sOff);

3233:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3234:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3235:        offset += dof;
3236:      }
3237:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3238: .ve

3240:   Level: developer

3242: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3243: @*/
3244: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3245: {
3246:   PetscSectionSym sym;
3247:   PetscErrorCode  ierr;

3252:   if (perms) *perms = NULL;
3253:   if (rots)  *rots  = NULL;
3254:   sym = section->sym;
3255:   if (sym && (perms || rots)) {
3256:     SymWorkLink link;

3258:     if (sym->workin) {
3259:       link        = sym->workin;
3260:       sym->workin = sym->workin->next;
3261:     } else {
3262:       PetscNewLog(sym,&link);
3263:     }
3264:     if (numPoints > link->numPoints) {
3265:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3266:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3267:       link->numPoints = numPoints;
3268:     }
3269:     link->next   = sym->workout;
3270:     sym->workout = link;
3271:     PetscArrayzero((PetscInt**)link->perms,numPoints);
3272:     PetscArrayzero((PetscInt**)link->rots,numPoints);
3273:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3274:     if (perms) *perms = link->perms;
3275:     if (rots)  *rots  = link->rots;
3276:   }
3277:   return(0);
3278: }

3280: /*@C
3281:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3283:   Not collective

3285:   Input Parameters:
3286: + section - the section
3287: . numPoints - the number of points
3288: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3289:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3290:     context, see DMPlexGetConeOrientation()).

3292:   Output Parameter:
3293: + perms - The permutations for the given orientations: set to NULL at conclusion
3294: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3296:   Level: developer

3298: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3299: @*/
3300: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3301: {
3302:   PetscSectionSym sym;

3306:   sym = section->sym;
3307:   if (sym && (perms || rots)) {
3308:     SymWorkLink *p,link;

3310:     for (p=&sym->workout; (link=*p); p=&link->next) {
3311:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3312:         *p          = link->next;
3313:         link->next  = sym->workin;
3314:         sym->workin = link;
3315:         if (perms) *perms = NULL;
3316:         if (rots)  *rots  = NULL;
3317:         return(0);
3318:       }
3319:     }
3320:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3321:   }
3322:   return(0);
3323: }

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

3328:   Not collective

3330:   Input Parameters:
3331: + section - the section
3332: . field - the field of the section
3333: . numPoints - the number of points
3334: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3335:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3336:     context, see DMPlexGetConeOrientation()).

3338:   Output Parameter:
3339: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3340: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3341:     identity).

3343:   Level: developer

3345: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3346: @*/
3347: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3348: {

3353:   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);
3354:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3355:   return(0);
3356: }

3358: /*@C
3359:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3361:   Not collective

3363:   Input Parameters:
3364: + section - the section
3365: . field - the field number
3366: . numPoints - the number of points
3367: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3368:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3369:     context, see DMPlexGetConeOrientation()).

3371:   Output Parameter:
3372: + perms - The permutations for the given orientations: set to NULL at conclusion
3373: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3375:   Level: developer

3377: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3378: @*/
3379: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3380: {

3385:   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);
3386:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3387:   return(0);
3388: }

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

3393:   Not collective

3395:   Input Parameter:
3396: . s - the global PetscSection

3398:   Output Parameters:
3399: . flg - the flag

3401:   Level: developer

3403: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3404: @*/
3405: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3406: {
3409:   *flg = s->useFieldOff;
3410:   return(0);
3411: }

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

3416:   Not collective

3418:   Input Parameters:
3419: + s   - the global PetscSection
3420: - flg - the flag

3422:   Level: developer

3424: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3425: @*/
3426: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3427: {
3430:   s->useFieldOff = flg;
3431:   return(0);
3432: }

3434: #define PetscSectionExpandPoints_Loop(TYPE) \
3435: { \
3436:   PetscInt i, n, o0, o1, size; \
3437:   TYPE *a0 = (TYPE*)origArray, *a1; \
3438:   PetscSectionGetStorageSize(s, &size); \
3439:   PetscMalloc1(size, &a1); \
3440:   for (i=0; i<npoints; i++) { \
3441:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3442:     PetscSectionGetOffset(s, i, &o1); \
3443:     PetscSectionGetDof(s, i, &n); \
3444:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3445:   } \
3446:   *newArray = (void*)a1; \
3447: }

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

3452:   Not collective

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

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

3464:   Level: developer

3466: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3467: @*/
3468: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3469: {
3470:   PetscSection        s;
3471:   const PetscInt      *points_;
3472:   PetscInt            i, n, npoints, pStart, pEnd;
3473:   PetscMPIInt         unitsize;
3474:   PetscErrorCode      ierr;

3482:   MPI_Type_size(dataType, &unitsize);
3483:   ISGetLocalSize(points, &npoints);
3484:   ISGetIndices(points, &points_);
3485:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3486:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3487:   PetscSectionSetChart(s, 0, npoints);
3488:   for (i=0; i<npoints; i++) {
3489:     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);
3490:     PetscSectionGetDof(origSection, points_[i], &n);
3491:     PetscSectionSetDof(s, i, n);
3492:   }
3493:   PetscSectionSetUp(s);
3494:   if (newArray) {
3495:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3496:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3497:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3498:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3499:   }
3500:   if (newSection) {
3501:     *newSection = s;
3502:   } else {
3503:     PetscSectionDestroy(&s);
3504:   }
3505:   return(0);
3506: }