Actual source code: section.c

petsc-master 2020-06-03
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)->compNames          = NULL;
 63:   (*s)->clObj              = NULL;
 64:   (*s)->clSection          = NULL;
 65:   (*s)->clPoints           = NULL;
 66:   (*s)->clSize             = 0;
 67:   (*s)->clPerm             = NULL;
 68:   (*s)->clInvPerm          = NULL;
 69:   return(0);
 70: }

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

 75:   Collective

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

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

 83:   Level: intermediate

 85: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
 86: @*/
 87: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 88: {
 89:   PetscSectionSym sym;
 90:   IS              perm;
 91:   PetscInt        numFields, f, c, pStart, pEnd, p;
 92:   PetscErrorCode  ierr;

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

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

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

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

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

164:   Collective

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

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

172:   Level: beginner

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

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

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

191:   Collective on PetscSection

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

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

199:   Level: intermediate

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

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

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

221:   Collective on PetscSection

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

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

230:   Level: intermediate

232:   Notes:
233:   Field names are disregarded.

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

250:   flg = PETSC_FALSE;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

326:   Not collective

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

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

334:   Level: intermediate

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

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

350:   Not collective

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

356:   Level: intermediate

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

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

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

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

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

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

393:   Not Collective

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

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

402:   Level: intermediate

404: .seealso: PetscSectionSetFieldName()
405: @*/
406: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
407: {
411:   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);
412:   *fieldName = s->fieldNames[field];
413:   return(0);
414: }

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

419:   Not Collective

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

426:   Level: intermediate

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

437:   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);
438:   PetscFree(s->fieldNames[field]);
439:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
440:   return(0);
441: }

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

446:   Not Collective

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

454:   Level: intermediate

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

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

472:   Not Collective

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

480:   Level: intermediate

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

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

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

501:   Not collective

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

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

510:   Level: intermediate

512: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
513: @*/
514: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
515: {
519:   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);
520:   *numComp = s->numFieldComponents[field];
521:   return(0);
522: }

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

527:   Not collective

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

534:   Level: intermediate

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

546:   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);
547:   if (s->compNames) {
548:     for (c = 0; c < s->numFieldComponents[field]; ++c) {
549:       PetscFree(s->compNames[field][c]);
550:     }
551:     PetscFree(s->compNames[field]);
552:   }

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

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

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

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

580:   Not collective

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

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

589:   Level: intermediate

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

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

605:   Not collective

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

612:   Level: intermediate

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

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

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

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

642:   Not collective

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

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

650:   Level: intermediate

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

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

665:   Not collective

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

671:   Level: intermediate

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

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

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

696:   Not collective

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

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

704:   Level: intermediate

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

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

720:   Not collective

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

726:   Not collective

728:   Level: intermediate

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

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

744:   Not collective

746:   Input Parameters:
747: + s - the PetscSection
748: - point - the point

750:   Output Parameter:
751: . numDof - the number of dof

753:   Level: intermediate

755: .seealso: PetscSectionSetDof(), PetscSectionCreate()
756: @*/
757: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
758: {
762:   if (PetscDefined(USE_DEBUG)) {
763:     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);
764:   }
765:   *numDof = s->atlasDof[point - s->pStart];
766:   return(0);
767: }

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

772:   Not collective

774:   Input Parameters:
775: + s - the PetscSection
776: . point - the point
777: - numDof - the number of dof

779:   Level: intermediate

781: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
782: @*/
783: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
784: {
787:   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);
788:   s->atlasDof[point - s->pStart] = numDof;
789:   return(0);
790: }

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

795:   Not collective

797:   Input Parameters:
798: + s - the PetscSection
799: . point - the point
800: - numDof - the number of additional dof

802:   Level: intermediate

804: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
805: @*/
806: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
807: {
810:   if (PetscDefined(USE_DEBUG)) {
811:     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);
812:   }
813:   s->atlasDof[point - s->pStart] += numDof;
814:   return(0);
815: }

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

820:   Not collective

822:   Input Parameters:
823: + s - the PetscSection
824: . point - the point
825: - field - the field

827:   Output Parameter:
828: . numDof - the number of dof

830:   Level: intermediate

832: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
833: @*/
834: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
835: {

840:   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);
841:   PetscSectionGetDof(s->field[field], point, numDof);
842:   return(0);
843: }

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

848:   Not collective

850:   Input Parameters:
851: + s - the PetscSection
852: . point - the point
853: . field - the field
854: - numDof - the number of dof

856:   Level: intermediate

858: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
859: @*/
860: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
861: {

866:   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);
867:   PetscSectionSetDof(s->field[field], point, numDof);
868:   return(0);
869: }

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

874:   Not collective

876:   Input Parameters:
877: + s - the PetscSection
878: . point - the point
879: . field - the field
880: - numDof - the number of dof

882:   Level: intermediate

884: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
885: @*/
886: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
887: {

892:   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);
893:   PetscSectionAddDof(s->field[field], point, numDof);
894:   return(0);
895: }

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

900:   Not collective

902:   Input Parameters:
903: + s - the PetscSection
904: - point - the point

906:   Output Parameter:
907: . numDof - the number of dof which are fixed by constraints

909:   Level: intermediate

911: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
912: @*/
913: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
914: {

920:   if (s->bc) {
921:     PetscSectionGetDof(s->bc, point, numDof);
922:   } else *numDof = 0;
923:   return(0);
924: }

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

929:   Not collective

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

936:   Level: intermediate

938: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
939: @*/
940: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
941: {

946:   if (numDof) {
947:     PetscSectionCheckConstraints_Static(s);
948:     PetscSectionSetDof(s->bc, point, numDof);
949:   }
950:   return(0);
951: }

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

956:   Not collective

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

963:   Level: intermediate

965: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
966: @*/
967: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
968: {

973:   if (numDof) {
974:     PetscSectionCheckConstraints_Static(s);
975:     PetscSectionAddDof(s->bc, point, numDof);
976:   }
977:   return(0);
978: }

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

983:   Not collective

985:   Input Parameters:
986: + s - the PetscSection
987: . point - the point
988: - field - the field

990:   Output Parameter:
991: . numDof - the number of dof which are fixed by constraints

993:   Level: intermediate

995: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
996: @*/
997: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
998: {

1004:   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);
1005:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
1006:   return(0);
1007: }

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

1012:   Not collective

1014:   Input Parameters:
1015: + s - the PetscSection
1016: . point - the point
1017: . field - the field
1018: - numDof - the number of dof which are fixed by constraints

1020:   Level: intermediate

1022: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1023: @*/
1024: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1025: {

1030:   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);
1031:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
1032:   return(0);
1033: }

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

1038:   Not collective

1040:   Input Parameters:
1041: + s - the PetscSection
1042: . point - the point
1043: . field - the field
1044: - numDof - the number of additional dof which are fixed by constraints

1046:   Level: intermediate

1048: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1049: @*/
1050: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1051: {

1056:   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);
1057:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
1058:   return(0);
1059: }

1061: /*@
1062:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1064:   Not collective

1066:   Input Parameter:
1067: . s - the PetscSection

1069:   Level: advanced

1071: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1072: @*/
1073: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1074: {

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

1082:     PetscSectionSetUp(s->bc);
1083:     PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1084:   }
1085:   return(0);
1086: }

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

1091:   Not collective

1093:   Input Parameter:
1094: . s - the PetscSection

1096:   Level: intermediate

1098: .seealso: PetscSectionCreate()
1099: @*/
1100: PetscErrorCode PetscSectionSetUp(PetscSection s)
1101: {
1102:   const PetscInt *pind   = NULL;
1103:   PetscInt        offset = 0, foff, p, f;
1104:   PetscErrorCode  ierr;

1108:   if (s->setup) return(0);
1109:   s->setup = PETSC_TRUE;
1110:   /* Set offsets and field offsets for all points */
1111:   /*   Assume that all fields have the same chart */
1112:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1113:   if (s->pointMajor) {
1114:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1115:       const PetscInt q = pind ? pind[p] : p;

1117:       /* Set point offset */
1118:       s->atlasOff[q] = offset;
1119:       offset        += s->atlasDof[q];
1120:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
1121:       /* Set field offset */
1122:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1123:         PetscSection sf = s->field[f];

1125:         sf->atlasOff[q] = foff;
1126:         foff += sf->atlasDof[q];
1127:       }
1128:     }
1129:   } else {
1130:     /* Set field offsets for all points */
1131:     for (f = 0; f < s->numFields; ++f) {
1132:       PetscSection sf = s->field[f];

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

1137:         sf->atlasOff[q] = offset;
1138:         offset += sf->atlasDof[q];
1139:       }
1140:     }
1141:     /* Disable point offsets since these are unused */
1142:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1143:       s->atlasOff[p] = -1;
1144:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
1145:     }
1146:   }
1147:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1148:   /* Setup BC sections */
1149:   PetscSectionSetUpBC(s);
1150:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1151:   return(0);
1152: }

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

1157:   Not collective

1159:   Input Parameters:
1160: . s - the PetscSection

1162:   Output Parameter:
1163: . maxDof - the maximum dof

1165:   Level: intermediate

1167: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1168: @*/
1169: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1170: {
1174:   *maxDof = s->maxDof;
1175:   return(0);
1176: }

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

1181:   Not collective

1183:   Input Parameter:
1184: . s - the PetscSection

1186:   Output Parameter:
1187: . size - the size of an array which can hold all the dofs

1189:   Level: intermediate

1191: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1192: @*/
1193: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1194: {
1195:   PetscInt p, n = 0;

1200:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1201:   *size = n;
1202:   return(0);
1203: }

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

1208:   Not collective

1210:   Input Parameter:
1211: . s - the PetscSection

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

1216:   Level: intermediate

1218: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1219: @*/
1220: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1221: {
1222:   PetscInt p, n = 0;

1227:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1228:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1229:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1230:   }
1231:   *size = n;
1232:   return(0);
1233: }

1235: /*@
1236:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1237:   the local section and an SF describing the section point overlap.

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

1245:   Output Parameter:
1246: . gsection - The PetscSection for the global field layout

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

1250:   Level: intermediate

1252: .seealso: PetscSectionCreate()
1253: @*/
1254: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1255: {
1256:   PetscSection    gs;
1257:   const PetscInt *pind = NULL;
1258:   PetscInt       *recv = NULL, *neg = NULL;
1259:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1260:   PetscErrorCode  ierr;

1268:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1269:   PetscSectionGetChart(s, &pStart, &pEnd);
1270:   PetscSectionSetChart(gs, pStart, pEnd);
1271:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1272:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1273:   /* Must allocate for all points visible to SF, which may be more than this section */
1274:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1275:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1276:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1277:     if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1278:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1279:     PetscArrayzero(neg,nroots);
1280:   }
1281:   /* Mark all local points with negative dof */
1282:   for (p = pStart; p < pEnd; ++p) {
1283:     PetscSectionGetDof(s, p, &dof);
1284:     PetscSectionSetDof(gs, p, dof);
1285:     PetscSectionGetConstraintDof(s, p, &cdof);
1286:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1287:     if (neg) neg[p] = -(dof+1);
1288:   }
1289:   PetscSectionSetUpBC(gs);
1290:   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]);}
1291:   if (nroots >= 0) {
1292:     PetscArrayzero(recv,nlocal);
1293:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1294:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1295:     for (p = pStart; p < pEnd; ++p) {
1296:       if (recv[p] < 0) {
1297:         gs->atlasDof[p-pStart] = recv[p];
1298:         PetscSectionGetDof(s, p, &dof);
1299:         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);
1300:       }
1301:     }
1302:   }
1303:   /* Calculate new sizes, get process offset, and calculate point offsets */
1304:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1305:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1306:     const PetscInt q = pind ? pind[p] : p;

1308:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1309:     gs->atlasOff[q] = off;
1310:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1311:   }
1312:   if (!localOffsets) {
1313:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1314:     globalOff -= off;
1315:   }
1316:   for (p = pStart, off = 0; p < pEnd; ++p) {
1317:     gs->atlasOff[p-pStart] += globalOff;
1318:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1319:   }
1320:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1321:   /* Put in negative offsets for ghost points */
1322:   if (nroots >= 0) {
1323:     PetscArrayzero(recv,nlocal);
1324:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1325:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1326:     for (p = pStart; p < pEnd; ++p) {
1327:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1328:     }
1329:   }
1330:   PetscFree2(neg,recv);
1331:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1332:   *gsection = gs;
1333:   return(0);
1334: }

1336: /*@
1337:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1338:   the local section and an SF describing the section point overlap.

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

1347:   Output Parameter:
1348:   . gsection - The PetscSection for the global field layout

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

1352:   Level: advanced

1354: .seealso: PetscSectionCreate()
1355: @*/
1356: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1357: {
1358:   const PetscInt *pind = NULL;
1359:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1360:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1361:   PetscErrorCode  ierr;

1367:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1368:   PetscSectionGetChart(s, &pStart, &pEnd);
1369:   PetscSectionSetChart(*gsection, pStart, pEnd);
1370:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1371:   if (nroots >= 0) {
1372:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1373:     PetscCalloc1(nroots, &neg);
1374:     if (nroots > pEnd-pStart) {
1375:       PetscCalloc1(nroots, &tmpOff);
1376:     } else {
1377:       tmpOff = &(*gsection)->atlasDof[-pStart];
1378:     }
1379:   }
1380:   /* Mark ghost points with negative dof */
1381:   for (p = pStart; p < pEnd; ++p) {
1382:     for (e = 0; e < numExcludes; ++e) {
1383:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1384:         PetscSectionSetDof(*gsection, p, 0);
1385:         break;
1386:       }
1387:     }
1388:     if (e < numExcludes) continue;
1389:     PetscSectionGetDof(s, p, &dof);
1390:     PetscSectionSetDof(*gsection, p, dof);
1391:     PetscSectionGetConstraintDof(s, p, &cdof);
1392:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1393:     if (neg) neg[p] = -(dof+1);
1394:   }
1395:   PetscSectionSetUpBC(*gsection);
1396:   if (nroots >= 0) {
1397:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1398:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1399:     if (nroots > pEnd - pStart) {
1400:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1401:     }
1402:   }
1403:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1404:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1405:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1406:     const PetscInt q = pind ? pind[p] : p;

1408:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1409:     (*gsection)->atlasOff[q] = off;
1410:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1411:   }
1412:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1413:   globalOff -= off;
1414:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1415:     (*gsection)->atlasOff[p] += globalOff;
1416:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1417:   }
1418:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1419:   /* Put in negative offsets for ghost points */
1420:   if (nroots >= 0) {
1421:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1422:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1423:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1424:     if (nroots > pEnd - pStart) {
1425:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1426:     }
1427:   }
1428:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1429:   PetscFree(neg);
1430:   return(0);
1431: }

1433: /*@
1434:   PetscSectionGetPointLayout - Get the PetscLayout associated with the section points

1436:   Collective on comm

1438:   Input Parameters:
1439: + comm - The MPI_Comm
1440: - s    - The PetscSection

1442:   Output Parameter:
1443: . layout - The point layout for the section

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

1447:   Level: advanced

1449: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1450: @*/
1451: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1452: {
1453:   PetscInt       pStart, pEnd, p, localSize = 0;

1457:   PetscSectionGetChart(s, &pStart, &pEnd);
1458:   for (p = pStart; p < pEnd; ++p) {
1459:     PetscInt dof;

1461:     PetscSectionGetDof(s, p, &dof);
1462:     if (dof > 0) ++localSize;
1463:   }
1464:   PetscLayoutCreate(comm, layout);
1465:   PetscLayoutSetLocalSize(*layout, localSize);
1466:   PetscLayoutSetBlockSize(*layout, 1);
1467:   PetscLayoutSetUp(*layout);
1468:   return(0);
1469: }

1471: /*@
1472:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

1474:   Collective on comm

1476:   Input Parameters:
1477: + comm - The MPI_Comm
1478: - s    - The PetscSection

1480:   Output Parameter:
1481: . layout - The dof layout for the section

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

1485:   Level: advanced

1487: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1488: @*/
1489: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1490: {
1491:   PetscInt       pStart, pEnd, p, localSize = 0;

1497:   PetscSectionGetChart(s, &pStart, &pEnd);
1498:   for (p = pStart; p < pEnd; ++p) {
1499:     PetscInt dof,cdof;

1501:     PetscSectionGetDof(s, p, &dof);
1502:     PetscSectionGetConstraintDof(s, p, &cdof);
1503:     if (dof-cdof > 0) localSize += dof-cdof;
1504:   }
1505:   PetscLayoutCreate(comm, layout);
1506:   PetscLayoutSetLocalSize(*layout, localSize);
1507:   PetscLayoutSetBlockSize(*layout, 1);
1508:   PetscLayoutSetUp(*layout);
1509:   return(0);
1510: }

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

1515:   Not collective

1517:   Input Parameters:
1518: + s - the PetscSection
1519: - point - the point

1521:   Output Parameter:
1522: . offset - the offset

1524:   Level: intermediate

1526: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1527: @*/
1528: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1529: {
1533:   if (PetscDefined(USE_DEBUG)) {
1534:     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);
1535:   }
1536:   *offset = s->atlasOff[point - s->pStart];
1537:   return(0);
1538: }

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

1543:   Not collective

1545:   Input Parameters:
1546: + s - the PetscSection
1547: . point - the point
1548: - offset - the offset

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

1552:   Level: intermediate

1554: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1555: @*/
1556: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1557: {
1560:   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);
1561:   s->atlasOff[point - s->pStart] = offset;
1562:   return(0);
1563: }

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

1568:   Not collective

1570:   Input Parameters:
1571: + s - the PetscSection
1572: . point - the point
1573: - field - the field

1575:   Output Parameter:
1576: . offset - the offset

1578:   Level: intermediate

1580: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1581: @*/
1582: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1583: {

1589:   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);
1590:   PetscSectionGetOffset(s->field[field], point, offset);
1591:   return(0);
1592: }

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

1597:   Not collective

1599:   Input Parameters:
1600: + s - the PetscSection
1601: . point - the point
1602: . field - the field
1603: - offset - the offset

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

1607:   Level: intermediate

1609: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1610: @*/
1611: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1612: {

1617:   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);
1618:   PetscSectionSetOffset(s->field[field], point, offset);
1619:   return(0);
1620: }

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

1625:   Not collective

1627:   Input Parameters:
1628: + s - the PetscSection
1629: . point - the point
1630: - field - the field

1632:   Output Parameter:
1633: . offset - the offset

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

1638:   Level: advanced

1640: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1641: @*/
1642: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1643: {
1644:   PetscInt       off, foff;

1650:   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);
1651:   PetscSectionGetOffset(s, point, &off);
1652:   PetscSectionGetOffset(s->field[field], point, &foff);
1653:   *offset = foff - off;
1654:   return(0);
1655: }

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

1660:   Not collective

1662:   Input Parameter:
1663: . s - the PetscSection

1665:   Output Parameters:
1666: + start - the minimum offset
1667: - end   - one more than the maximum offset

1669:   Level: intermediate

1671: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1672: @*/
1673: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1674: {
1675:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1680:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1681:   PetscSectionGetChart(s, &pStart, &pEnd);
1682:   for (p = 0; p < pEnd-pStart; ++p) {
1683:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1685:     if (off >= 0) {
1686:       os = PetscMin(os, off);
1687:       oe = PetscMax(oe, off+dof);
1688:     }
1689:   }
1690:   if (start) *start = os;
1691:   if (end)   *end   = oe;
1692:   return(0);
1693: }

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

1698:   Collective on s

1700:   Input Parameter:
1701: + s      - the PetscSection
1702: . len    - the number of subfields
1703: - fields - the subfield numbers

1705:   Output Parameter:
1706: . subs   - the subsection

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

1710:   Level: advanced

1712: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1713: @*/
1714: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1715: {
1716:   PetscInt       nF, f, c, pStart, pEnd, p, maxCdof = 0;

1720:   if (!len) return(0);
1724:   PetscSectionGetNumFields(s, &nF);
1725:   if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1726:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1727:   PetscSectionSetNumFields(*subs, len);
1728:   for (f = 0; f < len; ++f) {
1729:     const char *name   = NULL;
1730:     PetscInt   numComp = 0;

1732:     PetscSectionGetFieldName(s, fields[f], &name);
1733:     PetscSectionSetFieldName(*subs, f, name);
1734:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1735:     PetscSectionSetFieldComponents(*subs, f, numComp);
1736:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1737:       PetscSectionGetComponentName(s, fields[f], c, &name);
1738:       PetscSectionSetComponentName(*subs, f, c, name);
1739:     }
1740:   }
1741:   PetscSectionGetChart(s, &pStart, &pEnd);
1742:   PetscSectionSetChart(*subs, pStart, pEnd);
1743:   for (p = pStart; p < pEnd; ++p) {
1744:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1746:     for (f = 0; f < len; ++f) {
1747:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1748:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1749:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1750:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1751:       dof  += fdof;
1752:       cdof += cfdof;
1753:     }
1754:     PetscSectionSetDof(*subs, p, dof);
1755:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1756:     maxCdof = PetscMax(cdof, maxCdof);
1757:   }
1758:   PetscSectionSetUp(*subs);
1759:   if (maxCdof) {
1760:     PetscInt *indices;

1762:     PetscMalloc1(maxCdof, &indices);
1763:     for (p = pStart; p < pEnd; ++p) {
1764:       PetscInt cdof;

1766:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1767:       if (cdof) {
1768:         const PetscInt *oldIndices = NULL;
1769:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1774:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1775:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1776:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1777:           /* This can be sped up if we assume sorted fields */
1778:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1779:             PetscInt oldfdof = 0;
1780:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1781:             oldFoff += oldfdof;
1782:           }
1783:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1784:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1785:           numConst += cfdof;
1786:           fOff     += fdof;
1787:         }
1788:         PetscSectionSetConstraintIndices(*subs, p, indices);
1789:       }
1790:     }
1791:     PetscFree(indices);
1792:   }
1793:   return(0);
1794: }

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

1799:   Collective on s

1801:   Input Parameters:
1802: + s     - the input sections
1803: - len   - the number of input sections

1805:   Output Parameter:
1806: . supers - the supersection

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

1810:   Level: advanced

1812: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1813: @*/
1814: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1815: {
1816:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1820:   if (!len) return(0);
1821:   for (i = 0; i < len; ++i) {
1822:     PetscInt nf, pStarti, pEndi;

1824:     PetscSectionGetNumFields(s[i], &nf);
1825:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1826:     pStart = PetscMin(pStart, pStarti);
1827:     pEnd   = PetscMax(pEnd,   pEndi);
1828:     Nf += nf;
1829:   }
1830:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1831:   PetscSectionSetNumFields(*supers, Nf);
1832:   for (i = 0, f = 0; i < len; ++i) {
1833:     PetscInt nf, fi, ci;

1835:     PetscSectionGetNumFields(s[i], &nf);
1836:     for (fi = 0; fi < nf; ++fi, ++f) {
1837:       const char *name   = NULL;
1838:       PetscInt   numComp = 0;

1840:       PetscSectionGetFieldName(s[i], fi, &name);
1841:       PetscSectionSetFieldName(*supers, f, name);
1842:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1843:       PetscSectionSetFieldComponents(*supers, f, numComp);
1844:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1845:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1846:         PetscSectionSetComponentName(*supers, f, ci, name);
1847:       }
1848:     }
1849:   }
1850:   PetscSectionSetChart(*supers, pStart, pEnd);
1851:   for (p = pStart; p < pEnd; ++p) {
1852:     PetscInt dof = 0, cdof = 0;

1854:     for (i = 0, f = 0; i < len; ++i) {
1855:       PetscInt nf, fi, pStarti, pEndi;
1856:       PetscInt fdof = 0, cfdof = 0;

1858:       PetscSectionGetNumFields(s[i], &nf);
1859:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1860:       if ((p < pStarti) || (p >= pEndi)) continue;
1861:       for (fi = 0; fi < nf; ++fi, ++f) {
1862:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1863:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1864:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1865:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1866:         dof  += fdof;
1867:         cdof += cfdof;
1868:       }
1869:     }
1870:     PetscSectionSetDof(*supers, p, dof);
1871:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1872:     maxCdof = PetscMax(cdof, maxCdof);
1873:   }
1874:   PetscSectionSetUp(*supers);
1875:   if (maxCdof) {
1876:     PetscInt *indices;

1878:     PetscMalloc1(maxCdof, &indices);
1879:     for (p = pStart; p < pEnd; ++p) {
1880:       PetscInt cdof;

1882:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1883:       if (cdof) {
1884:         PetscInt dof, numConst = 0, fOff = 0;

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

1890:           PetscSectionGetNumFields(s[i], &nf);
1891:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1892:           if ((p < pStarti) || (p >= pEndi)) continue;
1893:           for (fi = 0; fi < nf; ++fi, ++f) {
1894:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1895:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1896:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1897:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1898:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1899:             numConst += cfdof;
1900:           }
1901:           PetscSectionGetDof(s[i], p, &dof);
1902:           fOff += dof;
1903:         }
1904:         PetscSectionSetConstraintIndices(*supers, p, indices);
1905:       }
1906:     }
1907:     PetscFree(indices);
1908:   }
1909:   return(0);
1910: }

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

1915:   Collective on s

1917:   Input Parameters:
1918: + s           - the PetscSection
1919: - subpointMap - a sorted list of points in the original mesh which are in the submesh

1921:   Output Parameter:
1922: . subs - the subsection

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

1926:   Level: advanced

1928: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1929: @*/
1930: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1931: {
1932:   const PetscInt *points = NULL, *indices = NULL;
1933:   PetscInt       numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;

1940:   PetscSectionGetNumFields(s, &numFields);
1941:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1942:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1943:   for (f = 0; f < numFields; ++f) {
1944:     const char *name   = NULL;
1945:     PetscInt   numComp = 0;

1947:     PetscSectionGetFieldName(s, f, &name);
1948:     PetscSectionSetFieldName(*subs, f, name);
1949:     PetscSectionGetFieldComponents(s, f, &numComp);
1950:     PetscSectionSetFieldComponents(*subs, f, numComp);
1951:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
1952:       PetscSectionGetComponentName(s, f, c, &name);
1953:       PetscSectionSetComponentName(*subs, f, c, name);
1954:     }
1955:   }
1956:   /* For right now, we do not try to squeeze the subchart */
1957:   if (subpointMap) {
1958:     ISGetSize(subpointMap, &numSubpoints);
1959:     ISGetIndices(subpointMap, &points);
1960:   }
1961:   PetscSectionGetChart(s, &pStart, &pEnd);
1962:   PetscSectionSetChart(*subs, 0, numSubpoints);
1963:   for (p = pStart; p < pEnd; ++p) {
1964:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1966:     PetscFindInt(p, numSubpoints, points, &subp);
1967:     if (subp < 0) continue;
1968:     for (f = 0; f < numFields; ++f) {
1969:       PetscSectionGetFieldDof(s, p, f, &fdof);
1970:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1971:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1972:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1973:     }
1974:     PetscSectionGetDof(s, p, &dof);
1975:     PetscSectionSetDof(*subs, subp, dof);
1976:     PetscSectionGetConstraintDof(s, p, &cdof);
1977:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1978:   }
1979:   PetscSectionSetUp(*subs);
1980:   /* Change offsets to original offsets */
1981:   for (p = pStart; p < pEnd; ++p) {
1982:     PetscInt off, foff = 0;

1984:     PetscFindInt(p, numSubpoints, points, &subp);
1985:     if (subp < 0) continue;
1986:     for (f = 0; f < numFields; ++f) {
1987:       PetscSectionGetFieldOffset(s, p, f, &foff);
1988:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1989:     }
1990:     PetscSectionGetOffset(s, p, &off);
1991:     PetscSectionSetOffset(*subs, subp, off);
1992:   }
1993:   /* Copy constraint indices */
1994:   for (subp = 0; subp < numSubpoints; ++subp) {
1995:     PetscInt cdof;

1997:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1998:     if (cdof) {
1999:       for (f = 0; f < numFields; ++f) {
2000:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
2001:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2002:       }
2003:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
2004:       PetscSectionSetConstraintIndices(*subs, subp, indices);
2005:     }
2006:   }
2007:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2008:   return(0);
2009: }

2011: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2012: {
2013:   PetscInt       p;
2014:   PetscMPIInt    rank;

2018:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2019:   PetscViewerASCIIPushSynchronized(viewer);
2020:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2021:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2022:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2023:       PetscInt b;

2025:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2026:       if (s->bcIndices) {
2027:         for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2028:           PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2029:         }
2030:       }
2031:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2032:     } else {
2033:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2034:     }
2035:   }
2036:   PetscViewerFlush(viewer);
2037:   PetscViewerASCIIPopSynchronized(viewer);
2038:   if (s->sym) {
2039:     PetscViewerASCIIPushTab(viewer);
2040:     PetscSectionSymView(s->sym,viewer);
2041:     PetscViewerASCIIPopTab(viewer);
2042:   }
2043:   return(0);
2044: }

2046: /*@C
2047:    PetscSectionViewFromOptions - View from Options

2049:    Collective on PetscSection

2051:    Input Parameters:
2052: +  A - the PetscSection object to view
2053: .  obj - Optional object
2054: -  name - command line option

2056:    Level: intermediate
2057: .seealso:  PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2058: @*/
2059: PetscErrorCode  PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2060: {

2065:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
2066:   return(0);
2067: }

2069: /*@C
2070:   PetscSectionView - Views a PetscSection

2072:   Collective on PetscSection

2074:   Input Parameters:
2075: + s - the PetscSection object to view
2076: - v - the viewer

2078:   Level: beginner

2080: .seealso PetscSectionCreate(), PetscSectionDestroy()
2081: @*/
2082: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2083: {
2084:   PetscBool      isascii;
2085:   PetscInt       f;

2090:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2092:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2093:   if (isascii) {
2094:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2095:     if (s->numFields) {
2096:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2097:       for (f = 0; f < s->numFields; ++f) {
2098:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
2099:         PetscSectionView_ASCII(s->field[f], viewer);
2100:       }
2101:     } else {
2102:       PetscSectionView_ASCII(s, viewer);
2103:     }
2104:   }
2105:   return(0);
2106: }

2108: /*@
2109:   PetscSectionReset - Frees all section data.

2111:   Not collective

2113:   Input Parameters:
2114: . s - the PetscSection

2116:   Level: beginner

2118: .seealso: PetscSection, PetscSectionCreate()
2119: @*/
2120: PetscErrorCode PetscSectionReset(PetscSection s)
2121: {
2122:   PetscInt       f, c;

2127:   for (f = 0; f < s->numFields; ++f) {
2128:     PetscSectionDestroy(&s->field[f]);
2129:     PetscFree(s->fieldNames[f]);
2130:     for (c = 0; c < s->numFieldComponents[f]; ++c)
2131:       PetscFree(s->compNames[f][c]);
2132:     PetscFree(s->compNames[f]);
2133:   }
2134:   PetscFree(s->numFieldComponents);
2135:   PetscFree(s->fieldNames);
2136:   PetscFree(s->compNames);
2137:   PetscFree(s->field);
2138:   PetscSectionDestroy(&s->bc);
2139:   PetscFree(s->bcIndices);
2140:   PetscFree2(s->atlasDof, s->atlasOff);
2141:   PetscSectionDestroy(&s->clSection);
2142:   ISDestroy(&s->clPoints);
2143:   ISDestroy(&s->perm);
2144:   PetscFree(s->clPerm);
2145:   PetscFree(s->clInvPerm);
2146:   PetscSectionSymDestroy(&s->sym);
2147:   PetscSectionDestroy(&s->clSection);
2148:   ISDestroy(&s->clPoints);

2150:   s->pStart    = -1;
2151:   s->pEnd      = -1;
2152:   s->maxDof    = 0;
2153:   s->setup     = PETSC_FALSE;
2154:   s->numFields = 0;
2155:   s->clObj     = NULL;
2156:   return(0);
2157: }

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

2162:   Not collective

2164:   Input Parameters:
2165: . s - the PetscSection

2167:   Level: beginner

2169: .seealso: PetscSection, PetscSectionCreate()
2170: @*/
2171: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2172: {

2176:   if (!*s) return(0);
2178:   if (--((PetscObject)(*s))->refct > 0) {
2179:     *s = NULL;
2180:     return(0);
2181:   }
2182:   PetscSectionReset(*s);
2183:   PetscHeaderDestroy(s);
2184:   return(0);
2185: }

2187: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2188: {
2189:   const PetscInt p = point - s->pStart;

2193:   *values = &baseArray[s->atlasOff[p]];
2194:   return(0);
2195: }

2197: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2198: {
2199:   PetscInt       *array;
2200:   const PetscInt p           = point - s->pStart;
2201:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2202:   PetscInt       cDim        = 0;

2207:   PetscSectionGetConstraintDof(s, p, &cDim);
2208:   array = &baseArray[s->atlasOff[p]];
2209:   if (!cDim) {
2210:     if (orientation >= 0) {
2211:       const PetscInt dim = s->atlasDof[p];
2212:       PetscInt       i;

2214:       if (mode == INSERT_VALUES) {
2215:         for (i = 0; i < dim; ++i) array[i] = values[i];
2216:       } else {
2217:         for (i = 0; i < dim; ++i) array[i] += values[i];
2218:       }
2219:     } else {
2220:       PetscInt offset = 0;
2221:       PetscInt j      = -1, field, i;

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

2226:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2227:         offset += dim;
2228:       }
2229:     }
2230:   } else {
2231:     if (orientation >= 0) {
2232:       const PetscInt dim  = s->atlasDof[p];
2233:       PetscInt       cInd = 0, i;
2234:       const PetscInt *cDof;

2236:       PetscSectionGetConstraintIndices(s, point, &cDof);
2237:       if (mode == INSERT_VALUES) {
2238:         for (i = 0; i < dim; ++i) {
2239:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2240:           array[i] = values[i];
2241:         }
2242:       } else {
2243:         for (i = 0; i < dim; ++i) {
2244:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2245:           array[i] += values[i];
2246:         }
2247:       }
2248:     } else {
2249:       const PetscInt *cDof;
2250:       PetscInt       offset  = 0;
2251:       PetscInt       cOffset = 0;
2252:       PetscInt       j       = 0, field;

2254:       PetscSectionGetConstraintIndices(s, point, &cDof);
2255:       for (field = 0; field < s->numFields; ++field) {
2256:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2257:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2258:         const PetscInt sDim = dim - tDim;
2259:         PetscInt       cInd = 0, i,k;

2261:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2262:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2263:           array[j] = values[k];
2264:         }
2265:         offset  += dim;
2266:         cOffset += dim - tDim;
2267:       }
2268:     }
2269:   }
2270:   return(0);
2271: }

2273: /*@C
2274:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2276:   Not collective

2278:   Input Parameter:
2279: . s - The PetscSection

2281:   Output Parameter:
2282: . hasConstraints - flag indicating that the section has constrained dofs

2284:   Level: intermediate

2286: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2287: @*/
2288: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2289: {
2293:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2294:   return(0);
2295: }

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

2300:   Not collective

2302:   Input Parameters:
2303: + s     - The PetscSection
2304: - point - The point

2306:   Output Parameter:
2307: . indices - The constrained dofs

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

2311:   Level: intermediate

2313: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2314: @*/
2315: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2316: {

2321:   if (s->bc) {
2322:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2323:   } else *indices = NULL;
2324:   return(0);
2325: }

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

2330:   Not collective

2332:   Input Parameters:
2333: + s     - The PetscSection
2334: . point - The point
2335: - indices - The constrained dofs

2337:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2339:   Level: intermediate

2341: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2342: @*/
2343: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2344: {

2349:   if (s->bc) {
2350:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2351:   }
2352:   return(0);
2353: }

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

2358:   Not collective

2360:   Input Parameters:
2361: + s     - The PetscSection
2362: . field  - The field number
2363: - point - The point

2365:   Output Parameter:
2366: . indices - The constrained dofs sorted in ascending order

2368:   Notes:
2369:   The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().

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

2374:   Level: intermediate

2376: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2377: @*/
2378: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2379: {

2384:   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);
2385:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2386:   return(0);
2387: }

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

2392:   Not collective

2394:   Input Parameters:
2395: + s       - The PetscSection
2396: . point   - The point
2397: . field   - The field number
2398: - indices - The constrained dofs

2400:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2402:   Level: intermediate

2404: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2405: @*/
2406: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2407: {

2412:   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);
2413:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2414:   return(0);
2415: }

2417: /*@
2418:   PetscSectionPermute - Reorder the section according to the input point permutation

2420:   Collective on PetscSection

2422:   Input Parameter:
2423: + section - The PetscSection object
2424: - perm - The point permutation, old point p becomes new point perm[p]

2426:   Output Parameter:
2427: . sectionNew - The permuted PetscSection

2429:   Level: intermediate

2431: .seealso: MatPermute()
2432: @*/
2433: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2434: {
2435:   PetscSection    s = section, sNew;
2436:   const PetscInt *perm;
2437:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;
2438:   PetscErrorCode  ierr;

2444:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2445:   PetscSectionGetNumFields(s, &numFields);
2446:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2447:   for (f = 0; f < numFields; ++f) {
2448:     const char *name;
2449:     PetscInt    numComp;

2451:     PetscSectionGetFieldName(s, f, &name);
2452:     PetscSectionSetFieldName(sNew, f, name);
2453:     PetscSectionGetFieldComponents(s, f, &numComp);
2454:     PetscSectionSetFieldComponents(sNew, f, numComp);
2455:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2456:       PetscSectionGetComponentName(s, f, c, &name);
2457:       PetscSectionSetComponentName(sNew, f, c, name);
2458:     }
2459:   }
2460:   ISGetLocalSize(permutation, &numPoints);
2461:   ISGetIndices(permutation, &perm);
2462:   PetscSectionGetChart(s, &pStart, &pEnd);
2463:   PetscSectionSetChart(sNew, pStart, pEnd);
2464:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2465:   for (p = pStart; p < pEnd; ++p) {
2466:     PetscInt dof, cdof;

2468:     PetscSectionGetDof(s, p, &dof);
2469:     PetscSectionSetDof(sNew, perm[p], dof);
2470:     PetscSectionGetConstraintDof(s, p, &cdof);
2471:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2472:     for (f = 0; f < numFields; ++f) {
2473:       PetscSectionGetFieldDof(s, p, f, &dof);
2474:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2475:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2476:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2477:     }
2478:   }
2479:   PetscSectionSetUp(sNew);
2480:   for (p = pStart; p < pEnd; ++p) {
2481:     const PetscInt *cind;
2482:     PetscInt        cdof;

2484:     PetscSectionGetConstraintDof(s, p, &cdof);
2485:     if (cdof) {
2486:       PetscSectionGetConstraintIndices(s, p, &cind);
2487:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2488:     }
2489:     for (f = 0; f < numFields; ++f) {
2490:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2491:       if (cdof) {
2492:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2493:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2494:       }
2495:     }
2496:   }
2497:   ISRestoreIndices(permutation, &perm);
2498:   *sectionNew = sNew;
2499:   return(0);
2500: }

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

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

2508:   Collective on sf

2510:   Input Parameters:
2511: + sf - The SF
2512: - rootSection - Section defined on root space

2514:   Output Parameters:
2515: + remoteOffsets - root offsets in leaf storage, or NULL
2516: - leafSection - Section defined on the leaf space

2518:   Level: advanced

2520: .seealso: PetscSFCreate()
2521: @*/
2522: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2523: {
2524:   PetscSF        embedSF;
2525:   const PetscInt *indices;
2526:   IS             selected;
2527:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
2528:   PetscBool      *sub, hasc;

2532:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2533:   PetscSectionGetNumFields(rootSection, &numFields);
2534:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2535:   PetscMalloc1(numFields+2, &sub);
2536:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2537:   for (f = 0; f < numFields; ++f) {
2538:     PetscSectionSym sym;
2539:     const char      *name   = NULL;
2540:     PetscInt        numComp = 0;

2542:     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2543:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2544:     PetscSectionGetFieldName(rootSection, f, &name);
2545:     PetscSectionGetFieldSym(rootSection, f, &sym);
2546:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2547:     PetscSectionSetFieldName(leafSection, f, name);
2548:     PetscSectionSetFieldSym(leafSection, f, sym);
2549:     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
2550:       PetscSectionGetComponentName(rootSection, f, c, &name);
2551:       PetscSectionSetComponentName(leafSection, f, c, name);
2552:     }
2553:   }
2554:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2555:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2556:   rpEnd = PetscMin(rpEnd,nroots);
2557:   rpEnd = PetscMax(rpStart,rpEnd);
2558:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2559:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2560:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2561:   if (sub[0]) {
2562:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2563:     ISGetIndices(selected, &indices);
2564:     PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2565:     ISRestoreIndices(selected, &indices);
2566:     ISDestroy(&selected);
2567:   } else {
2568:     PetscObjectReference((PetscObject)sf);
2569:     embedSF = sf;
2570:   }
2571:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2572:   lpEnd++;

2574:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

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

2580:   /* Could fuse these at the cost of copies and extra allocation */
2581:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2582:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2583:   if (sub[1]) {
2584:     PetscSectionCheckConstraints_Static(rootSection);
2585:     PetscSectionCheckConstraints_Static(leafSection);
2586:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2587:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2588:   }
2589:   for (f = 0; f < numFields; ++f) {
2590:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2591:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2592:     if (sub[2+f]) {
2593:       PetscSectionCheckConstraints_Static(rootSection->field[f]);
2594:       PetscSectionCheckConstraints_Static(leafSection->field[f]);
2595:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2596:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2597:     }
2598:   }
2599:   if (remoteOffsets) {
2600:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2601:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2602:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2603:   }
2604:   PetscSectionSetUp(leafSection);
2605:   if (hasc) { /* need to communicate bcIndices */
2606:     PetscSF  bcSF;
2607:     PetscInt *rOffBc;

2609:     PetscMalloc1(lpEnd - lpStart, &rOffBc);
2610:     if (sub[1]) {
2611:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2612:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2613:       PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2614:       PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2615:       PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2616:       PetscSFDestroy(&bcSF);
2617:     }
2618:     for (f = 0; f < numFields; ++f) {
2619:       if (sub[2+f]) {
2620:         PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2621:         PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2622:         PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2623:         PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2624:         PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2625:         PetscSFDestroy(&bcSF);
2626:       }
2627:     }
2628:     PetscFree(rOffBc);
2629:   }
2630:   PetscSFDestroy(&embedSF);
2631:   PetscFree(sub);
2632:   PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2633:   return(0);
2634: }

2636: /*@C
2637:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

2639:   Collective on sf

2641:   Input Parameters:
2642: + sf          - The SF
2643: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2644: - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)

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

2649:   Level: developer

2651: .seealso: PetscSFCreate()
2652: @*/
2653: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2654: {
2655:   PetscSF         embedSF;
2656:   const PetscInt *indices;
2657:   IS              selected;
2658:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2659:   PetscErrorCode  ierr;

2662:   *remoteOffsets = NULL;
2663:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2664:   if (numRoots < 0) return(0);
2665:   PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2666:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2667:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2668:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2669:   ISGetIndices(selected, &indices);
2670:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2671:   ISRestoreIndices(selected, &indices);
2672:   ISDestroy(&selected);
2673:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2674:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2675:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2676:   PetscSFDestroy(&embedSF);
2677:   PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2678:   return(0);
2679: }

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

2684:   Collective on sf

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

2692:   Output Parameters:
2693: - sectionSF - The new SF

2695:   Note: Either rootSection or remoteOffsets can be specified

2697:   Level: advanced

2699: .seealso: PetscSFCreate()
2700: @*/
2701: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2702: {
2703:   MPI_Comm          comm;
2704:   const PetscInt    *localPoints;
2705:   const PetscSFNode *remotePoints;
2706:   PetscInt          lpStart, lpEnd;
2707:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2708:   PetscInt          *localIndices;
2709:   PetscSFNode       *remoteIndices;
2710:   PetscInt          i, ind;
2711:   PetscErrorCode    ierr;

2719:   PetscObjectGetComm((PetscObject)sf,&comm);
2720:   PetscSFCreate(comm, sectionSF);
2721:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2722:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2723:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2724:   if (numRoots < 0) return(0);
2725:   PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2726:   for (i = 0; i < numPoints; ++i) {
2727:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2728:     PetscInt dof;

2730:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2731:       PetscSectionGetDof(leafSection, localPoint, &dof);
2732:       numIndices += dof;
2733:     }
2734:   }
2735:   PetscMalloc1(numIndices, &localIndices);
2736:   PetscMalloc1(numIndices, &remoteIndices);
2737:   /* Create new index graph */
2738:   for (i = 0, ind = 0; i < numPoints; ++i) {
2739:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2740:     PetscInt rank       = remotePoints[i].rank;

2742:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2743:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2744:       PetscInt loff, dof, d;

2746:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2747:       PetscSectionGetDof(leafSection, localPoint, &dof);
2748:       for (d = 0; d < dof; ++d, ++ind) {
2749:         localIndices[ind]        = loff+d;
2750:         remoteIndices[ind].rank  = rank;
2751:         remoteIndices[ind].index = remoteOffset+d;
2752:       }
2753:     }
2754:   }
2755:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2756:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2757:   PetscSFSetUp(*sectionSF);
2758:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2759:   return(0);
2760: }

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

2765:   Collective on section

2767:   Input Parameters:
2768: + section   - The PetscSection
2769: . obj       - A PetscObject which serves as the key for this index
2770: . clSection - Section giving the size of the closure of each point
2771: - clPoints  - IS giving the points in each closure

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

2775:   Level: advanced

2777: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2778: @*/
2779: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2780: {

2787:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2788:   section->clObj     = obj;
2789:   PetscObjectReference((PetscObject)clSection);
2790:   PetscObjectReference((PetscObject)clPoints);
2791:   PetscSectionDestroy(&section->clSection);
2792:   ISDestroy(&section->clPoints);
2793:   section->clSection = clSection;
2794:   section->clPoints  = clPoints;
2795:   return(0);
2796: }

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

2801:   Collective on section

2803:   Input Parameters:
2804: + section   - The PetscSection
2805: - obj       - A PetscObject which serves as the key for this index

2807:   Output Parameters:
2808: + clSection - Section giving the size of the closure of each point
2809: - clPoints  - IS giving the points in each closure

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

2813:   Level: advanced

2815: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2816: @*/
2817: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2818: {
2820:   if (section->clObj == obj) {
2821:     if (clSection) *clSection = section->clSection;
2822:     if (clPoints)  *clPoints  = section->clPoints;
2823:   } else {
2824:     if (clSection) *clSection = NULL;
2825:     if (clPoints)  *clPoints  = NULL;
2826:   }
2827:   return(0);
2828: }

2830: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2831: {
2832:   PetscInt       i;

2836:   if (section->clObj != obj) {
2837:     PetscSectionDestroy(&section->clSection);
2838:     ISDestroy(&section->clPoints);
2839:   }
2840:   section->clObj  = obj;
2841:   PetscFree(section->clPerm);
2842:   PetscFree(section->clInvPerm);
2843:   section->clSize = clSize;
2844:   if (mode == PETSC_COPY_VALUES) {
2845:     PetscMalloc1(clSize, &section->clPerm);
2846:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2847:     PetscArraycpy(section->clPerm, clPerm, clSize);
2848:   } else if (mode == PETSC_OWN_POINTER) {
2849:     section->clPerm = clPerm;
2850:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2851:   PetscMalloc1(clSize, &section->clInvPerm);
2852:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2853:   return(0);
2854: }

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

2859:   Not Collective

2861:   Input Parameters:
2862: + section - The PetscSection
2863: . obj     - A PetscObject which serves as the key for this index
2864: - perm    - Permutation of the cell dof closure

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

2869:   Level: intermediate

2871: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2872: @*/
2873: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2874: {
2875:   const PetscInt *clPerm = NULL;
2876:   PetscInt        clSize = 0;
2877:   PetscErrorCode  ierr;

2880:   if (perm) {
2881:     ISGetLocalSize(perm, &clSize);
2882:     ISGetIndices(perm, &clPerm);
2883:   }
2884:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2885:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2886:   return(0);
2887: }

2889: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2890: {
2892:   if (section->clObj == obj) {
2893:     if (size) *size = section->clSize;
2894:     if (perm) *perm = section->clPerm;
2895:   } else {
2896:     if (size) *size = 0;
2897:     if (perm) *perm = NULL;
2898:   }
2899:   return(0);
2900: }

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

2905:   Not collective

2907:   Input Parameters:
2908: + section   - The PetscSection
2909: - obj       - A PetscObject which serves as the key for this index

2911:   Output Parameter:
2912: . perm - The dof closure permutation

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

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

2919:   Level: intermediate

2921: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2922: @*/
2923: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2924: {
2925:   const PetscInt *clPerm;
2926:   PetscInt        clSize;
2927:   PetscErrorCode  ierr;

2930:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2931:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2932:   return(0);
2933: }

2935: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2936: {
2938:   if (section->clObj == obj) {
2939:     if (size) *size = section->clSize;
2940:     if (perm) *perm = section->clInvPerm;
2941:   } else {
2942:     if (size) *size = 0;
2943:     if (perm) *perm = NULL;
2944:   }
2945:   return(0);
2946: }

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

2951:   Not collective

2953:   Input Parameters:
2954: + section   - The PetscSection
2955: - obj       - A PetscObject which serves as the key for this index

2957:   Output Parameters:
2958: + size - The dof closure size
2959: - perm - The dof closure permutation

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

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

2966:   Level: intermediate

2968: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2969: @*/
2970: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2971: {
2972:   const PetscInt *clPerm;
2973:   PetscInt        clSize;
2974:   PetscErrorCode  ierr;

2977:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2978:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2979:   return(0);
2980: }

2982: /*@
2983:   PetscSectionGetField - Get the subsection associated with a single field

2985:   Input Parameters:
2986: + s     - The PetscSection
2987: - field - The field number

2989:   Output Parameter:
2990: . subs  - The subsection for the given field

2992:   Level: intermediate

2994: .seealso: PetscSectionSetNumFields()
2995: @*/
2996: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2997: {
3001:   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);
3002:   *subs = s->field[field];
3003:   return(0);
3004: }

3006: PetscClassId      PETSC_SECTION_SYM_CLASSID;
3007: PetscFunctionList PetscSectionSymList = NULL;

3009: /*@
3010:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

3012:   Collective

3014:   Input Parameter:
3015: . comm - the MPI communicator

3017:   Output Parameter:
3018: . sym - pointer to the new set of symmetries

3020:   Level: developer

3022: .seealso: PetscSectionSym, PetscSectionSymDestroy()
3023: @*/
3024: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3025: {

3030:   ISInitializePackage();
3031:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
3032:   return(0);
3033: }

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

3038:   Collective on PetscSectionSym

3040:   Input Parameters:
3041: + sym    - The section symmetry object
3042: - method - The name of the section symmetry type

3044:   Level: developer

3046: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
3047: @*/
3048: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3049: {
3050:   PetscErrorCode (*r)(PetscSectionSym);
3051:   PetscBool      match;

3056:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
3057:   if (match) return(0);

3059:   PetscFunctionListFind(PetscSectionSymList,method,&r);
3060:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3061:   if (sym->ops->destroy) {
3062:     (*sym->ops->destroy)(sym);
3063:     sym->ops->destroy = NULL;
3064:   }
3065:   (*r)(sym);
3066:   PetscObjectChangeTypeName((PetscObject)sym,method);
3067:   return(0);
3068: }


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

3074:   Not Collective

3076:   Input Parameter:
3077: . sym  - The section symmetry

3079:   Output Parameter:
3080: . type - The index set type name

3082:   Level: developer

3084: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
3085: @*/
3086: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3087: {
3091:   *type = ((PetscObject)sym)->type_name;
3092:   return(0);
3093: }

3095: /*@C
3096:   PetscSectionSymRegister - Adds a new section symmetry implementation

3098:   Not Collective

3100:   Input Parameters:
3101: + name        - The name of a new user-defined creation routine
3102: - create_func - The creation routine itself

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

3107:   Level: developer

3109: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3110: @*/
3111: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3112: {

3116:   ISInitializePackage();
3117:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3118:   return(0);
3119: }

3121: /*@
3122:    PetscSectionSymDestroy - Destroys a section symmetry.

3124:    Collective on PetscSectionSym

3126:    Input Parameters:
3127: .  sym - the section symmetry

3129:    Level: developer

3131: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3132: @*/
3133: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3134: {
3135:   SymWorkLink    link,next;

3139:   if (!*sym) return(0);
3141:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
3142:   if ((*sym)->ops->destroy) {
3143:     (*(*sym)->ops->destroy)(*sym);
3144:   }
3145:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3146:   for (link=(*sym)->workin; link; link=next) {
3147:     next = link->next;
3148:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3149:     PetscFree(link);
3150:   }
3151:   (*sym)->workin = NULL;
3152:   PetscHeaderDestroy(sym);
3153:   return(0);
3154: }

3156: /*@C
3157:    PetscSectionSymView - Displays a section symmetry

3159:    Collective on PetscSectionSym

3161:    Input Parameters:
3162: +  sym - the index set
3163: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

3165:    Level: developer

3167: .seealso: PetscViewerASCIIOpen()
3168: @*/
3169: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3170: {

3175:   if (!viewer) {
3176:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3177:   }
3180:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3181:   if (sym->ops->view) {
3182:     (*sym->ops->view)(sym,viewer);
3183:   }
3184:   return(0);
3185: }

3187: /*@
3188:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3190:   Collective on PetscSection

3192:   Input Parameters:
3193: + section - the section describing data layout
3194: - sym - the symmetry describing the affect of orientation on the access of the data

3196:   Level: developer

3198: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3199: @*/
3200: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3201: {

3206:   PetscSectionSymDestroy(&(section->sym));
3207:   if (sym) {
3210:     PetscObjectReference((PetscObject) sym);
3211:   }
3212:   section->sym = sym;
3213:   return(0);
3214: }

3216: /*@
3217:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3219:   Not collective

3221:   Input Parameters:
3222: . section - the section describing data layout

3224:   Output Parameters:
3225: . sym - the symmetry describing the affect of orientation on the access of the data

3227:   Level: developer

3229: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3230: @*/
3231: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3232: {
3235:   *sym = section->sym;
3236:   return(0);
3237: }

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

3242:   Collective on PetscSection

3244:   Input Parameters:
3245: + section - the section describing data layout
3246: . field - the field number
3247: - sym - the symmetry describing the affect of orientation on the access of the data

3249:   Level: developer

3251: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3252: @*/
3253: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3254: {

3259:   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);
3260:   PetscSectionSetSym(section->field[field],sym);
3261:   return(0);
3262: }

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

3267:   Collective on PetscSection

3269:   Input Parameters:
3270: + section - the section describing data layout
3271: - field - the field number

3273:   Output Parameters:
3274: . sym - the symmetry describing the affect of orientation on the access of the data

3276:   Level: developer

3278: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3279: @*/
3280: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3281: {
3284:   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);
3285:   *sym = section->field[field]->sym;
3286:   return(0);
3287: }

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

3292:   Not collective

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

3301:   Output Parameter:
3302: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3303: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3304:     identity).

3306:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3307: .vb
3308:      const PetscInt    **perms;
3309:      const PetscScalar **rots;
3310:      PetscInt            lOffset;

3312:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3313:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3314:        PetscInt           point = points[2*i], dof, sOffset;
3315:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3316:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3318:        PetscSectionGetDof(section,point,&dof);
3319:        PetscSectionGetOffset(section,point,&sOffset);

3321:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3322:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3323:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3324:        lOffset += dof;
3325:      }
3326:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3327: .ve

3329:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3330: .vb
3331:      const PetscInt    **perms;
3332:      const PetscScalar **rots;
3333:      PetscInt            lOffset;

3335:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3336:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3337:        PetscInt           point = points[2*i], dof, sOffset;
3338:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3339:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3341:        PetscSectionGetDof(section,point,&dof);
3342:        PetscSectionGetOffset(section,point,&sOff);

3344:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3345:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3346:        offset += dof;
3347:      }
3348:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3349: .ve

3351:   Level: developer

3353: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3354: @*/
3355: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3356: {
3357:   PetscSectionSym sym;
3358:   PetscErrorCode  ierr;

3363:   if (perms) *perms = NULL;
3364:   if (rots)  *rots  = NULL;
3365:   sym = section->sym;
3366:   if (sym && (perms || rots)) {
3367:     SymWorkLink link;

3369:     if (sym->workin) {
3370:       link        = sym->workin;
3371:       sym->workin = sym->workin->next;
3372:     } else {
3373:       PetscNewLog(sym,&link);
3374:     }
3375:     if (numPoints > link->numPoints) {
3376:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3377:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3378:       link->numPoints = numPoints;
3379:     }
3380:     link->next   = sym->workout;
3381:     sym->workout = link;
3382:     PetscArrayzero((PetscInt**)link->perms,numPoints);
3383:     PetscArrayzero((PetscInt**)link->rots,numPoints);
3384:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3385:     if (perms) *perms = link->perms;
3386:     if (rots)  *rots  = link->rots;
3387:   }
3388:   return(0);
3389: }

3391: /*@C
3392:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3394:   Not collective

3396:   Input Parameters:
3397: + section - the section
3398: . numPoints - the number of points
3399: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3400:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3401:     context, see DMPlexGetConeOrientation()).

3403:   Output Parameter:
3404: + perms - The permutations for the given orientations: set to NULL at conclusion
3405: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3407:   Level: developer

3409: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3410: @*/
3411: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3412: {
3413:   PetscSectionSym sym;

3417:   sym = section->sym;
3418:   if (sym && (perms || rots)) {
3419:     SymWorkLink *p,link;

3421:     for (p=&sym->workout; (link=*p); p=&link->next) {
3422:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3423:         *p          = link->next;
3424:         link->next  = sym->workin;
3425:         sym->workin = link;
3426:         if (perms) *perms = NULL;
3427:         if (rots)  *rots  = NULL;
3428:         return(0);
3429:       }
3430:     }
3431:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3432:   }
3433:   return(0);
3434: }

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

3439:   Not collective

3441:   Input Parameters:
3442: + section - the section
3443: . field - the field of the section
3444: . numPoints - the number of points
3445: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3446:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3447:     context, see DMPlexGetConeOrientation()).

3449:   Output Parameter:
3450: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3451: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3452:     identity).

3454:   Level: developer

3456: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3457: @*/
3458: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3459: {

3464:   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);
3465:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3466:   return(0);
3467: }

3469: /*@C
3470:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3472:   Not collective

3474:   Input Parameters:
3475: + section - the section
3476: . field - the field number
3477: . numPoints - the number of points
3478: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3479:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3480:     context, see DMPlexGetConeOrientation()).

3482:   Output Parameter:
3483: + perms - The permutations for the given orientations: set to NULL at conclusion
3484: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3486:   Level: developer

3488: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3489: @*/
3490: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3491: {

3496:   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);
3497:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3498:   return(0);
3499: }

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

3504:   Not collective

3506:   Input Parameter:
3507: . s - the global PetscSection

3509:   Output Parameters:
3510: . flg - the flag

3512:   Level: developer

3514: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3515: @*/
3516: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3517: {
3520:   *flg = s->useFieldOff;
3521:   return(0);
3522: }

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

3527:   Not collective

3529:   Input Parameters:
3530: + s   - the global PetscSection
3531: - flg - the flag

3533:   Level: developer

3535: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3536: @*/
3537: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3538: {
3541:   s->useFieldOff = flg;
3542:   return(0);
3543: }

3545: #define PetscSectionExpandPoints_Loop(TYPE) \
3546: { \
3547:   PetscInt i, n, o0, o1, size; \
3548:   TYPE *a0 = (TYPE*)origArray, *a1; \
3549:   PetscSectionGetStorageSize(s, &size); \
3550:   PetscMalloc1(size, &a1); \
3551:   for (i=0; i<npoints; i++) { \
3552:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3553:     PetscSectionGetOffset(s, i, &o1); \
3554:     PetscSectionGetDof(s, i, &n); \
3555:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3556:   } \
3557:   *newArray = (void*)a1; \
3558: }

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

3563:   Not collective

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

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

3575:   Level: developer

3577: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3578: @*/
3579: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3580: {
3581:   PetscSection        s;
3582:   const PetscInt      *points_;
3583:   PetscInt            i, n, npoints, pStart, pEnd;
3584:   PetscMPIInt         unitsize;
3585:   PetscErrorCode      ierr;

3593:   MPI_Type_size(dataType, &unitsize);
3594:   ISGetLocalSize(points, &npoints);
3595:   ISGetIndices(points, &points_);
3596:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3597:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3598:   PetscSectionSetChart(s, 0, npoints);
3599:   for (i=0; i<npoints; i++) {
3600:     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);
3601:     PetscSectionGetDof(origSection, points_[i], &n);
3602:     PetscSectionSetDof(s, i, n);
3603:   }
3604:   PetscSectionSetUp(s);
3605:   if (newArray) {
3606:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3607:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3608:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3609:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3610:   }
3611:   if (newSection) {
3612:     *newSection = s;
3613:   } else {
3614:     PetscSectionDestroy(&s);
3615:   }
3616:   return(0);
3617: }