Actual source code: vsectionis.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

  5:  #include <petsc/private/isimpl.h>
  6:  #include <petscsf.h>
  7:  #include <petscviewer.h>

  9: PetscClassId PETSC_SECTION_CLASSID;

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

 14:   Collective on MPI_Comm

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

 20:   Level: developer

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

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

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

 43:   ISInitializePackage();

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

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

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

 73:   Collective on MPI_Comm

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

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

 81:   Level: developer

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

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

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

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

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

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

157:   Collective on MPI_Comm

159:   Input Parameter:
160: . section - the PetscSection

162:   Output Parameter:
163: . newSection - the copy

165:   Level: developer

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

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

181: /*@
182:   PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database

184:   Collective on PetscSection

186:   Input Parameter:
187: . section - the PetscSection

189:   Options Database:
190: . -petscsection_point_major the dof order

192:   Level: intermediate

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

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

211: /*@
212:   PetscSectionCompare - Compares two sections

214:   Collective on PetscSection

216:   Input Parameterss:
217: + s1 - the first PetscSection
218: - s2 - the second PetscSection

220:   Output Parameter:
221: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise 

223:   Level: developer

225:   Notes:
226:   Field names are disregarded.

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

243:   flg = PETSC_FALSE;

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

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

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

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

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

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

275:     PetscSectionGetConstraintIndices(s1, p, &idx1);
276:     PetscSectionGetConstraintIndices(s2, p, &idx2);
277:     PetscMemcmp(idx1, idx2, ncdof*sizeof(PetscInt), congruent);
278:     if (!(*congruent)) goto not_congruent;
279:   }

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

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

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

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

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

303:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
304:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
305:       PetscMemcmp(idx1, idx2, nfcdof*sizeof(PetscInt), congruent);
306:       if (!(*congruent)) goto not_congruent;
307:     }
308:   }

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

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

319:   Not collective

321:   Input Parameter:
322: . s - the PetscSection

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

327:   Level: intermediate

329: .seealso: PetscSectionSetNumFields()
330: @*/
331: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
332: {
336:   *numFields = s->numFields;
337:   return(0);
338: }

340: /*@
341:   PetscSectionSetNumFields - Sets the number of fields.

343:   Not collective

345:   Input Parameters:
346: + s - the PetscSection
347: - numFields - the number of fields

349:   Level: intermediate

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

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

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

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

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

379: /*@C
380:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

382:   Not Collective

384:   Input Parameters:
385: + s     - the PetscSection
386: - field - the field number

388:   Output Parameter:
389: . fieldName - the field name

391:   Level: developer

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

405: /*@C
406:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

408:   Not Collective

410:   Input Parameters:
411: + s     - the PetscSection
412: . field - the field number
413: - fieldName - the field name

415:   Level: developer

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

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

432: /*@
433:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

435:   Not collective

437:   Input Parameters:
438: + s - the PetscSection
439: - field - the field number

441:   Output Parameter:
442: . numComp - the number of field components

444:   Level: intermediate

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

458: /*@
459:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

461:   Not collective

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

468:   Level: intermediate

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

481: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
482: {

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

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

496:   Not collective

498:   Input Parameter:
499: . s - the PetscSection

501:   Output Parameters:
502: + pStart - the first point
503: - pEnd - one past the last point

505:   Level: intermediate

507: .seealso: PetscSectionSetChart(), PetscSectionCreate()
508: @*/
509: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
510: {
513:   if (pStart) *pStart = s->pStart;
514:   if (pEnd)   *pEnd   = s->pEnd;
515:   return(0);
516: }

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

521:   Not collective

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

528:   Level: intermediate

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

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

545:   s->pStart = pStart;
546:   s->pEnd   = pEnd;
547:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
548:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
549:   for (f = 0; f < s->numFields; ++f) {
550:     PetscSectionSetChart(s->field[f], pStart, pEnd);
551:   }
552:   return(0);
553: }

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

558:   Not collective

560:   Input Parameter:
561: . s - the PetscSection

563:   Output Parameters:
564: . perm - The permutation as an IS

566:   Level: intermediate

568: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
569: @*/
570: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
571: {
575:   return(0);
576: }

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

581:   Not collective

583:   Input Parameters:
584: + s - the PetscSection
585: - perm - the permutation of points

587:   Level: intermediate

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

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

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

612:   Not collective

614:   Input Parameter:
615: . s - the PetscSection

617:   Output Parameter:
618: . pm - the flag for point major ordering

620:   Level: intermediate

622: .seealso: PetscSectionSetPointMajor()
623: @*/
624: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
625: {
629:   *pm = s->pointMajor;
630:   return(0);
631: }

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

636:   Not collective

638:   Input Parameters:
639: + s  - the PetscSection
640: - pm - the flag for point major ordering

642:   Not collective

644:   Level: intermediate

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

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

660:   Not collective

662:   Input Parameters:
663: + s - the PetscSection
664: - point - the point

666:   Output Parameter:
667: . numDof - the number of dof

669:   Level: intermediate

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

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

688:   Not collective

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

695:   Level: intermediate

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

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

711:   Not collective

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

718:   Level: intermediate

720: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
721: @*/
722: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
723: {
726:   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);
727:   s->atlasDof[point - s->pStart] += numDof;
728:   return(0);
729: }

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

734:   Not collective

736:   Input Parameters:
737: + s - the PetscSection
738: . point - the point
739: - field - the field

741:   Output Parameter:
742: . numDof - the number of dof

744:   Level: intermediate

746: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
747: @*/
748: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
749: {

754:   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);
755:   PetscSectionGetDof(s->field[field], point, numDof);
756:   return(0);
757: }

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

762:   Not collective

764:   Input Parameters:
765: + s - the PetscSection
766: . point - the point
767: . field - the field
768: - numDof - the number of dof

770:   Level: intermediate

772: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
773: @*/
774: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
775: {

780:   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);
781:   PetscSectionSetDof(s->field[field], point, numDof);
782:   return(0);
783: }

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

788:   Not collective

790:   Input Parameters:
791: + s - the PetscSection
792: . point - the point
793: . field - the field
794: - numDof - the number of dof

796:   Level: intermediate

798: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
799: @*/
800: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
801: {

806:   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);
807:   PetscSectionAddDof(s->field[field], point, numDof);
808:   return(0);
809: }

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

814:   Not collective

816:   Input Parameters:
817: + s - the PetscSection
818: - point - the point

820:   Output Parameter:
821: . numDof - the number of dof which are fixed by constraints

823:   Level: intermediate

825: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
826: @*/
827: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
828: {

834:   if (s->bc) {
835:     PetscSectionGetDof(s->bc, point, numDof);
836:   } else *numDof = 0;
837:   return(0);
838: }

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

843:   Not collective

845:   Input Parameters:
846: + s - the PetscSection
847: . point - the point
848: - numDof - the number of dof which are fixed by constraints

850:   Level: intermediate

852: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
853: @*/
854: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
855: {

860:   if (numDof) {
861:     PetscSectionCheckConstraints_Static(s);
862:     PetscSectionSetDof(s->bc, point, numDof);
863:   }
864:   return(0);
865: }

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

870:   Not collective

872:   Input Parameters:
873: + s - the PetscSection
874: . point - the point
875: - numDof - the number of additional dof which are fixed by constraints

877:   Level: intermediate

879: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
880: @*/
881: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
882: {

887:   if (numDof) {
888:     PetscSectionCheckConstraints_Static(s);
889:     PetscSectionAddDof(s->bc, point, numDof);
890:   }
891:   return(0);
892: }

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

897:   Not collective

899:   Input Parameters:
900: + s - the PetscSection
901: . point - the point
902: - field - the field

904:   Output Parameter:
905: . numDof - the number of dof which are fixed by constraints

907:   Level: intermediate

909: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
910: @*/
911: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
912: {

918:   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);
919:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
920:   return(0);
921: }

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

926:   Not collective

928:   Input Parameters:
929: + s - the PetscSection
930: . point - the point
931: . field - the field
932: - numDof - the number of dof which are fixed by constraints

934:   Level: intermediate

936: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
937: @*/
938: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
939: {

944:   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);
945:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
946:   return(0);
947: }

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

952:   Not collective

954:   Input Parameters:
955: + s - the PetscSection
956: . point - the point
957: . field - the field
958: - numDof - the number of additional dof which are fixed by constraints

960:   Level: intermediate

962: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
963: @*/
964: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
965: {

970:   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);
971:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
972:   return(0);
973: }

975: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
976: {

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

984:     PetscSectionSetUp(s->bc);
985:     PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
986:   }
987:   return(0);
988: }

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

993:   Not collective

995:   Input Parameter:
996: . s - the PetscSection

998:   Level: intermediate

1000: .seealso: PetscSectionCreate()
1001: @*/
1002: PetscErrorCode PetscSectionSetUp(PetscSection s)
1003: {
1004:   const PetscInt *pind   = NULL;
1005:   PetscInt        offset = 0, foff, p, f;
1006:   PetscErrorCode  ierr;

1010:   if (s->setup) return(0);
1011:   s->setup = PETSC_TRUE;
1012:   /* Set offsets and field offsets for all points */
1013:   /*   Assume that all fields have the same chart */
1014:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1015:   if (s->pointMajor) {
1016:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1017:       const PetscInt q = pind ? pind[p] : p;

1019:       /* Set point offset */
1020:       s->atlasOff[q] = offset;
1021:       offset        += s->atlasDof[q];
1022:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
1023:       /* Set field offset */
1024:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1025:         PetscSection sf = s->field[f];

1027:         sf->atlasOff[q] = foff;
1028:         foff += sf->atlasDof[q];
1029:       }
1030:     }
1031:   } else {
1032:     /* Set field offsets for all points */
1033:     for (f = 0; f < s->numFields; ++f) {
1034:       PetscSection sf = s->field[f];

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

1039:         sf->atlasOff[q] = offset;
1040:         offset += sf->atlasDof[q];
1041:       }
1042:     }
1043:     /* Disable point offsets since these are unused */
1044:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1045:       s->atlasOff[p] = -1;
1046:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
1047:     }
1048:   }
1049:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1050:   /* Setup BC sections */
1051:   PetscSectionSetUpBC(s);
1052:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1053:   return(0);
1054: }

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

1059:   Not collective

1061:   Input Parameters:
1062: . s - the PetscSection

1064:   Output Parameter:
1065: . maxDof - the maximum dof

1067:   Level: intermediate

1069: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1070: @*/
1071: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1072: {
1076:   *maxDof = s->maxDof;
1077:   return(0);
1078: }

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

1083:   Not collective

1085:   Input Parameters:
1086: + s - the PetscSection
1087: - size - the allocated size

1089:   Output Parameter:
1090: . size - the size of an array which can hold all the dofs

1092:   Level: intermediate

1094: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1095: @*/
1096: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1097: {
1098:   PetscInt p, n = 0;

1103:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1104:   *size = n;
1105:   return(0);
1106: }

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

1111:   Not collective

1113:   Input Parameters:
1114: + s - the PetscSection
1115: - point - the point

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

1120:   Level: intermediate

1122: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1123: @*/
1124: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1125: {
1126:   PetscInt p, n = 0;

1131:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1132:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1133:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1134:   }
1135:   *size = n;
1136:   return(0);
1137: }

1139: /*@
1140:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1141:   the local section and an SF describing the section point overlap.

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

1149:   Output Parameter:
1150:   . gsection - The PetscSection for the global field layout

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

1154:   Level: developer

1156: .seealso: PetscSectionCreate()
1157: @*/
1158: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1159: {
1160:   PetscSection    gs;
1161:   const PetscInt *pind = NULL;
1162:   PetscInt       *recv = NULL, *neg = NULL;
1163:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
1164:   PetscErrorCode  ierr;

1172:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1173:   PetscSectionGetChart(s, &pStart, &pEnd);
1174:   PetscSectionSetChart(gs, pStart, pEnd);
1175:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1176:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1177:   /* Must allocate for all points visible to SF, which may be more than this section */
1178:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1179:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1180:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1181:     PetscMemzero(neg,nroots*sizeof(PetscInt));
1182:   }
1183:   /* Mark all local points with negative dof */
1184:   for (p = pStart; p < pEnd; ++p) {
1185:     PetscSectionGetDof(s, p, &dof);
1186:     PetscSectionSetDof(gs, p, dof);
1187:     PetscSectionGetConstraintDof(s, p, &cdof);
1188:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1189:     if (neg) neg[p] = -(dof+1);
1190:   }
1191:   PetscSectionSetUpBC(gs);
1192:   if (gs->bcIndices) {PetscMemcpy(gs->bcIndices, s->bcIndices, sizeof(PetscInt) * (gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]));}
1193:   if (nroots >= 0) {
1194:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1195:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1196:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1197:     for (p = pStart; p < pEnd; ++p) {
1198:       if (recv[p] < 0) {
1199:         gs->atlasDof[p-pStart] = recv[p];
1200:         PetscSectionGetDof(s, p, &dof);
1201:         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);
1202:       }
1203:     }
1204:   }
1205:   /* Calculate new sizes, get process offset, and calculate point offsets */
1206:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1207:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1208:     const PetscInt q = pind ? pind[p] : p;

1210:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1211:     gs->atlasOff[q] = off;
1212:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1213:   }
1214:   if (!localOffsets) {
1215:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1216:     globalOff -= off;
1217:   }
1218:   for (p = pStart, off = 0; p < pEnd; ++p) {
1219:     gs->atlasOff[p-pStart] += globalOff;
1220:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1221:   }
1222:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1223:   /* Put in negative offsets for ghost points */
1224:   if (nroots >= 0) {
1225:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1226:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1227:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1228:     for (p = pStart; p < pEnd; ++p) {
1229:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1230:     }
1231:   }
1232:   PetscFree2(neg,recv);
1233:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1234:   *gsection = gs;
1235:   return(0);
1236: }

1238: /*@
1239:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1240:   the local section and an SF describing the section point overlap.

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

1249:   Output Parameter:
1250:   . gsection - The PetscSection for the global field layout

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

1254:   Level: developer

1256: .seealso: PetscSectionCreate()
1257: @*/
1258: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1259: {
1260:   const PetscInt *pind = NULL;
1261:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1262:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1263:   PetscErrorCode  ierr;

1269:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1270:   PetscSectionGetChart(s, &pStart, &pEnd);
1271:   PetscSectionSetChart(*gsection, pStart, pEnd);
1272:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1273:   if (nroots >= 0) {
1274:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1275:     PetscCalloc1(nroots, &neg);
1276:     if (nroots > pEnd-pStart) {
1277:       PetscCalloc1(nroots, &tmpOff);
1278:     } else {
1279:       tmpOff = &(*gsection)->atlasDof[-pStart];
1280:     }
1281:   }
1282:   /* Mark ghost points with negative dof */
1283:   for (p = pStart; p < pEnd; ++p) {
1284:     for (e = 0; e < numExcludes; ++e) {
1285:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1286:         PetscSectionSetDof(*gsection, p, 0);
1287:         break;
1288:       }
1289:     }
1290:     if (e < numExcludes) continue;
1291:     PetscSectionGetDof(s, p, &dof);
1292:     PetscSectionSetDof(*gsection, p, dof);
1293:     PetscSectionGetConstraintDof(s, p, &cdof);
1294:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1295:     if (neg) neg[p] = -(dof+1);
1296:   }
1297:   PetscSectionSetUpBC(*gsection);
1298:   if (nroots >= 0) {
1299:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1300:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1301:     if (nroots > pEnd - pStart) {
1302:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1303:     }
1304:   }
1305:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1306:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1307:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1308:     const PetscInt q = pind ? pind[p] : p;

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

1335: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1336: {
1337:   PetscInt       pStart, pEnd, p, localSize = 0;

1341:   PetscSectionGetChart(s, &pStart, &pEnd);
1342:   for (p = pStart; p < pEnd; ++p) {
1343:     PetscInt dof;

1345:     PetscSectionGetDof(s, p, &dof);
1346:     if (dof > 0) ++localSize;
1347:   }
1348:   PetscLayoutCreate(comm, layout);
1349:   PetscLayoutSetLocalSize(*layout, localSize);
1350:   PetscLayoutSetBlockSize(*layout, 1);
1351:   PetscLayoutSetUp(*layout);
1352:   return(0);
1353: }

1355: /*@
1356:   PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.

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

1362:   Output Parameter:
1363: . layout - The layout for the section

1365:   Level: developer

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

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

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

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

1395:   Not collective

1397:   Input Parameters:
1398: + s - the PetscSection
1399: - point - the point

1401:   Output Parameter:
1402: . offset - the offset

1404:   Level: intermediate

1406: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1407: @*/
1408: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1409: {
1413: #if defined(PETSC_USE_DEBUG)
1414:   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);
1415: #endif
1416:   *offset = s->atlasOff[point - s->pStart];
1417:   return(0);
1418: }

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

1423:   Not collective

1425:   Input Parameters:
1426: + s - the PetscSection
1427: . point - the point
1428: - offset - the offset

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

1432:   Level: intermediate

1434: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1435: @*/
1436: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1437: {
1440:   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);
1441:   s->atlasOff[point - s->pStart] = offset;
1442:   return(0);
1443: }

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

1448:   Not collective

1450:   Input Parameters:
1451: + s - the PetscSection
1452: . point - the point
1453: - field - the field

1455:   Output Parameter:
1456: . offset - the offset

1458:   Level: intermediate

1460: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1461: @*/
1462: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1463: {

1469:   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);
1470:   PetscSectionGetOffset(s->field[field], point, offset);
1471:   return(0);
1472: }

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

1477:   Not collective

1479:   Input Parameters:
1480: + s - the PetscSection
1481: . point - the point
1482: . field - the field
1483: - offset - the offset

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

1487:   Level: intermediate

1489: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1490: @*/
1491: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1492: {

1497:   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);
1498:   PetscSectionSetOffset(s->field[field], point, offset);
1499:   return(0);
1500: }

1502: /* This gives the offset on a point of the field, ignoring constraints */
1503: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1504: {
1505:   PetscInt       off, foff;

1511:   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);
1512:   PetscSectionGetOffset(s, point, &off);
1513:   PetscSectionGetOffset(s->field[field], point, &foff);
1514:   *offset = foff - off;
1515:   return(0);
1516: }

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

1521:   Not collective

1523:   Input Parameter:
1524: . s - the PetscSection

1526:   Output Parameters:
1527: + start - the minimum offset
1528: - end   - one more than the maximum offset

1530:   Level: intermediate

1532: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1533: @*/
1534: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1535: {
1536:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1541:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1542:   PetscSectionGetChart(s, &pStart, &pEnd);
1543:   for (p = 0; p < pEnd-pStart; ++p) {
1544:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1546:     if (off >= 0) {
1547:       os = PetscMin(os, off);
1548:       oe = PetscMax(oe, off+dof);
1549:     }
1550:   }
1551:   if (start) *start = os;
1552:   if (end)   *end   = oe;
1553:   return(0);
1554: }

1556: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, const PetscInt fields[], PetscSection *subs)
1557: {
1558:   PetscInt       nF, f, pStart, pEnd, p, maxCdof = 0;

1562:   if (!numFields) return(0);
1566:   PetscSectionGetNumFields(s, &nF);
1567:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", numFields, nF);
1568:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1569:   PetscSectionSetNumFields(*subs, numFields);
1570:   for (f = 0; f < numFields; ++f) {
1571:     const char *name   = NULL;
1572:     PetscInt   numComp = 0;

1574:     PetscSectionGetFieldName(s, fields[f], &name);
1575:     PetscSectionSetFieldName(*subs, f, name);
1576:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1577:     PetscSectionSetFieldComponents(*subs, f, numComp);
1578:   }
1579:   PetscSectionGetChart(s, &pStart, &pEnd);
1580:   PetscSectionSetChart(*subs, pStart, pEnd);
1581:   for (p = pStart; p < pEnd; ++p) {
1582:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1584:     for (f = 0; f < numFields; ++f) {
1585:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1586:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1587:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1588:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1589:       dof  += fdof;
1590:       cdof += cfdof;
1591:     }
1592:     PetscSectionSetDof(*subs, p, dof);
1593:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1594:     maxCdof = PetscMax(cdof, maxCdof);
1595:   }
1596:   PetscSectionSetUp(*subs);
1597:   if (maxCdof) {
1598:     PetscInt *indices;

1600:     PetscMalloc1(maxCdof, &indices);
1601:     for (p = pStart; p < pEnd; ++p) {
1602:       PetscInt cdof;

1604:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1605:       if (cdof) {
1606:         const PetscInt *oldIndices = NULL;
1607:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1609:         for (f = 0; f < numFields; ++f) {
1610:           PetscInt oldFoff = 0, oldf;

1612:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1613:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1614:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1615:           /* This can be sped up if we assume sorted fields */
1616:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1617:             PetscInt oldfdof = 0;
1618:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1619:             oldFoff += oldfdof;
1620:           }
1621:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1622:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1623:           numConst += cfdof;
1624:           fOff     += fdof;
1625:         }
1626:         PetscSectionSetConstraintIndices(*subs, p, indices);
1627:       }
1628:     }
1629:     PetscFree(indices);
1630:   }
1631:   return(0);
1632: }

1634: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1635: {
1636:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1640:   if (!len) return(0);
1641:   for (i = 0; i < len; ++i) {
1642:     PetscInt nf, pStarti, pEndi;

1644:     PetscSectionGetNumFields(s[i], &nf);
1645:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1646:     pStart = PetscMin(pStart, pStarti);
1647:     pEnd   = PetscMax(pEnd,   pEndi);
1648:     Nf += nf;
1649:   }
1650:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1651:   PetscSectionSetNumFields(*supers, Nf);
1652:   for (i = 0, f = 0; i < len; ++i) {
1653:     PetscInt nf, fi;

1655:     PetscSectionGetNumFields(s[i], &nf);
1656:     for (fi = 0; fi < nf; ++fi, ++f) {
1657:       const char *name   = NULL;
1658:       PetscInt   numComp = 0;

1660:       PetscSectionGetFieldName(s[i], fi, &name);
1661:       PetscSectionSetFieldName(*supers, f, name);
1662:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1663:       PetscSectionSetFieldComponents(*supers, f, numComp);
1664:     }
1665:   }
1666:   PetscSectionSetChart(*supers, pStart, pEnd);
1667:   for (p = pStart; p < pEnd; ++p) {
1668:     PetscInt dof = 0, cdof = 0;

1670:     for (i = 0, f = 0; i < len; ++i) {
1671:       PetscInt nf, fi, pStarti, pEndi;
1672:       PetscInt fdof = 0, cfdof = 0;

1674:       PetscSectionGetNumFields(s[i], &nf);
1675:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1676:       if ((p < pStarti) || (p >= pEndi)) continue;
1677:       for (fi = 0; fi < nf; ++fi, ++f) {
1678:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1679:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1680:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1681:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1682:         dof  += fdof;
1683:         cdof += cfdof;
1684:       }
1685:     }
1686:     PetscSectionSetDof(*supers, p, dof);
1687:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1688:     maxCdof = PetscMax(cdof, maxCdof);
1689:   }
1690:   PetscSectionSetUp(*supers);
1691:   if (maxCdof) {
1692:     PetscInt *indices;

1694:     PetscMalloc1(maxCdof, &indices);
1695:     for (p = pStart; p < pEnd; ++p) {
1696:       PetscInt cdof;

1698:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1699:       if (cdof) {
1700:         PetscInt dof, numConst = 0, fOff = 0;

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

1706:           PetscSectionGetNumFields(s[i], &nf);
1707:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1708:           if ((p < pStarti) || (p >= pEndi)) continue;
1709:           for (fi = 0; fi < nf; ++fi, ++f) {
1710:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1711:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1712:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1713:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1714:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1715:             numConst += cfdof;
1716:           }
1717:           PetscSectionGetDof(s[i], p, &dof);
1718:           fOff += dof;
1719:         }
1720:         PetscSectionSetConstraintIndices(*supers, p, indices);
1721:       }
1722:     }
1723:     PetscFree(indices);
1724:   }
1725:   return(0);
1726: }

1728: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1729: {
1730:   const PetscInt *points = NULL, *indices = NULL;
1731:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1738:   PetscSectionGetNumFields(s, &numFields);
1739:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1740:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1741:   for (f = 0; f < numFields; ++f) {
1742:     const char *name   = NULL;
1743:     PetscInt   numComp = 0;

1745:     PetscSectionGetFieldName(s, f, &name);
1746:     PetscSectionSetFieldName(*subs, f, name);
1747:     PetscSectionGetFieldComponents(s, f, &numComp);
1748:     PetscSectionSetFieldComponents(*subs, f, numComp);
1749:   }
1750:   /* For right now, we do not try to squeeze the subchart */
1751:   if (subpointMap) {
1752:     ISGetSize(subpointMap, &numSubpoints);
1753:     ISGetIndices(subpointMap, &points);
1754:   }
1755:   PetscSectionGetChart(s, &pStart, &pEnd);
1756:   PetscSectionSetChart(*subs, 0, numSubpoints);
1757:   for (p = pStart; p < pEnd; ++p) {
1758:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1760:     PetscFindInt(p, numSubpoints, points, &subp);
1761:     if (subp < 0) continue;
1762:     for (f = 0; f < numFields; ++f) {
1763:       PetscSectionGetFieldDof(s, p, f, &fdof);
1764:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1765:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1766:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1767:     }
1768:     PetscSectionGetDof(s, p, &dof);
1769:     PetscSectionSetDof(*subs, subp, dof);
1770:     PetscSectionGetConstraintDof(s, p, &cdof);
1771:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1772:   }
1773:   PetscSectionSetUp(*subs);
1774:   /* Change offsets to original offsets */
1775:   for (p = pStart; p < pEnd; ++p) {
1776:     PetscInt off, foff = 0;

1778:     PetscFindInt(p, numSubpoints, points, &subp);
1779:     if (subp < 0) continue;
1780:     for (f = 0; f < numFields; ++f) {
1781:       PetscSectionGetFieldOffset(s, p, f, &foff);
1782:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1783:     }
1784:     PetscSectionGetOffset(s, p, &off);
1785:     PetscSectionSetOffset(*subs, subp, off);
1786:   }
1787:   /* Copy constraint indices */
1788:   for (subp = 0; subp < numSubpoints; ++subp) {
1789:     PetscInt cdof;

1791:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1792:     if (cdof) {
1793:       for (f = 0; f < numFields; ++f) {
1794:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1795:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1796:       }
1797:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1798:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1799:     }
1800:   }
1801:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1802:   return(0);
1803: }

1805: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1806: {
1807:   PetscInt       p;
1808:   PetscMPIInt    rank;

1812:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1813:   PetscViewerASCIIPushSynchronized(viewer);
1814:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1815:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1816:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1817:       PetscInt b;

1819:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1820:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1821:         PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
1822:       }
1823:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1824:     } else {
1825:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1826:     }
1827:   }
1828:   PetscViewerFlush(viewer);
1829:   PetscViewerASCIIPopSynchronized(viewer);
1830:   if (s->sym) {
1831:     PetscViewerASCIIPushTab(viewer);
1832:     PetscSectionSymView(s->sym,viewer);
1833:     PetscViewerASCIIPopTab(viewer);
1834:   }
1835:   return(0);
1836: }

1838: /*@C
1839:   PetscSectionView - Views a PetscSection

1841:   Collective on PetscSection

1843:   Input Parameters:
1844: + s - the PetscSection object to view
1845: - v - the viewer

1847:   Level: developer

1849: .seealso PetscSectionCreate(), PetscSectionDestroy()
1850: @*/
1851: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1852: {
1853:   PetscBool      isascii;
1854:   PetscInt       f;

1859:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1861:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1862:   if (isascii) {
1863:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1864:     if (s->numFields) {
1865:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1866:       for (f = 0; f < s->numFields; ++f) {
1867:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1868:         PetscSectionView_ASCII(s->field[f], viewer);
1869:       }
1870:     } else {
1871:       PetscSectionView_ASCII(s, viewer);
1872:     }
1873:   }
1874:   return(0);
1875: }

1877: /*@
1878:   PetscSectionReset - Frees all section data.

1880:   Not collective

1882:   Input Parameters:
1883: . s - the PetscSection

1885:   Level: developer

1887: .seealso: PetscSection, PetscSectionCreate()
1888: @*/
1889: PetscErrorCode PetscSectionReset(PetscSection s)
1890: {
1891:   PetscInt       f;

1896:   PetscFree(s->numFieldComponents);
1897:   for (f = 0; f < s->numFields; ++f) {
1898:     PetscSectionDestroy(&s->field[f]);
1899:     PetscFree(s->fieldNames[f]);
1900:   }
1901:   PetscFree(s->fieldNames);
1902:   PetscFree(s->field);
1903:   PetscSectionDestroy(&s->bc);
1904:   PetscFree(s->bcIndices);
1905:   PetscFree2(s->atlasDof, s->atlasOff);
1906:   PetscSectionDestroy(&s->clSection);
1907:   ISDestroy(&s->clPoints);
1908:   ISDestroy(&s->perm);
1909:   PetscFree(s->clPerm);
1910:   PetscFree(s->clInvPerm);
1911:   PetscSectionSymDestroy(&s->sym);

1913:   s->pStart    = -1;
1914:   s->pEnd      = -1;
1915:   s->maxDof    = 0;
1916:   s->setup     = PETSC_FALSE;
1917:   s->numFields = 0;
1918:   s->clObj     = NULL;
1919:   return(0);
1920: }

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

1925:   Not collective

1927:   Input Parameters:
1928: . s - the PetscSection

1930:   Level: developer

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

1935: .seealso: PetscSection, PetscSectionCreate()
1936: @*/
1937: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1938: {

1942:   if (!*s) return(0);
1944:   if (--((PetscObject)(*s))->refct > 0) {
1945:     *s = NULL;
1946:     return(0);
1947:   }
1948:   PetscSectionReset(*s);
1949:   PetscHeaderDestroy(s);
1950:   return(0);
1951: }

1953: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1954: {
1955:   const PetscInt p = point - s->pStart;

1959:   *values = &baseArray[s->atlasOff[p]];
1960:   return(0);
1961: }

1963: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1964: {
1965:   PetscInt       *array;
1966:   const PetscInt p           = point - s->pStart;
1967:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1968:   PetscInt       cDim        = 0;

1973:   PetscSectionGetConstraintDof(s, p, &cDim);
1974:   array = &baseArray[s->atlasOff[p]];
1975:   if (!cDim) {
1976:     if (orientation >= 0) {
1977:       const PetscInt dim = s->atlasDof[p];
1978:       PetscInt       i;

1980:       if (mode == INSERT_VALUES) {
1981:         for (i = 0; i < dim; ++i) array[i] = values[i];
1982:       } else {
1983:         for (i = 0; i < dim; ++i) array[i] += values[i];
1984:       }
1985:     } else {
1986:       PetscInt offset = 0;
1987:       PetscInt j      = -1, field, i;

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

1992:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1993:         offset += dim;
1994:       }
1995:     }
1996:   } else {
1997:     if (orientation >= 0) {
1998:       const PetscInt dim  = s->atlasDof[p];
1999:       PetscInt       cInd = 0, i;
2000:       const PetscInt *cDof;

2002:       PetscSectionGetConstraintIndices(s, point, &cDof);
2003:       if (mode == INSERT_VALUES) {
2004:         for (i = 0; i < dim; ++i) {
2005:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2006:           array[i] = values[i];
2007:         }
2008:       } else {
2009:         for (i = 0; i < dim; ++i) {
2010:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2011:           array[i] += values[i];
2012:         }
2013:       }
2014:     } else {
2015:       const PetscInt *cDof;
2016:       PetscInt       offset  = 0;
2017:       PetscInt       cOffset = 0;
2018:       PetscInt       j       = 0, field;

2020:       PetscSectionGetConstraintIndices(s, point, &cDof);
2021:       for (field = 0; field < s->numFields; ++field) {
2022:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2023:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2024:         const PetscInt sDim = dim - tDim;
2025:         PetscInt       cInd = 0, i,k;

2027:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2028:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2029:           array[j] = values[k];
2030:         }
2031:         offset  += dim;
2032:         cOffset += dim - tDim;
2033:       }
2034:     }
2035:   }
2036:   return(0);
2037: }

2039: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2040: {
2044:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2045:   return(0);
2046: }

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

2051:   Input Parameters:
2052: + s     - The PetscSection
2053: - point - The point

2055:   Output Parameter:
2056: . indices - The constrained dofs

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

2060:   Level: advanced

2062: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2063: @*/
2064: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2065: {

2070:   if (s->bc) {
2071:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2072:   } else *indices = NULL;
2073:   return(0);
2074: }

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

2079:   Input Parameters:
2080: + s     - The PetscSection
2081: . point - The point
2082: - indices - The constrained dofs

2084:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2086:   Level: advanced

2088: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2089: @*/
2090: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2091: {

2096:   if (s->bc) {
2097:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2098:   }
2099:   return(0);
2100: }

2102: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2103: {

2108:   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);
2109:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2110:   return(0);
2111: }

2113: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2114: {

2119:   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);
2120:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2121:   return(0);
2122: }

2124: /*@
2125:   PetscSectionPermute - Reorder the section according to the input point permutation

2127:   Collective on PetscSection

2129:   Input Parameter:
2130: + section - The PetscSection object
2131: - perm - The point permutation, old point p becomes new point perm[p]

2133:   Output Parameter:
2134: . sectionNew - The permuted PetscSection

2136:   Level: intermediate

2138: .keywords: mesh
2139: .seealso: MatPermute()
2140: @*/
2141: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2142: {
2143:   PetscSection    s = section, sNew;
2144:   const PetscInt *perm;
2145:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
2146:   PetscErrorCode  ierr;

2152:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2153:   PetscSectionGetNumFields(s, &numFields);
2154:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2155:   for (f = 0; f < numFields; ++f) {
2156:     const char *name;
2157:     PetscInt    numComp;

2159:     PetscSectionGetFieldName(s, f, &name);
2160:     PetscSectionSetFieldName(sNew, f, name);
2161:     PetscSectionGetFieldComponents(s, f, &numComp);
2162:     PetscSectionSetFieldComponents(sNew, f, numComp);
2163:   }
2164:   ISGetLocalSize(permutation, &numPoints);
2165:   ISGetIndices(permutation, &perm);
2166:   PetscSectionGetChart(s, &pStart, &pEnd);
2167:   PetscSectionSetChart(sNew, pStart, pEnd);
2168:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2169:   for (p = pStart; p < pEnd; ++p) {
2170:     PetscInt dof, cdof;

2172:     PetscSectionGetDof(s, p, &dof);
2173:     PetscSectionSetDof(sNew, perm[p], dof);
2174:     PetscSectionGetConstraintDof(s, p, &cdof);
2175:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2176:     for (f = 0; f < numFields; ++f) {
2177:       PetscSectionGetFieldDof(s, p, f, &dof);
2178:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2179:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2180:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2181:     }
2182:   }
2183:   PetscSectionSetUp(sNew);
2184:   for (p = pStart; p < pEnd; ++p) {
2185:     const PetscInt *cind;
2186:     PetscInt        cdof;

2188:     PetscSectionGetConstraintDof(s, p, &cdof);
2189:     if (cdof) {
2190:       PetscSectionGetConstraintIndices(s, p, &cind);
2191:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2192:     }
2193:     for (f = 0; f < numFields; ++f) {
2194:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2195:       if (cdof) {
2196:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2197:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2198:       }
2199:     }
2200:   }
2201:   ISRestoreIndices(permutation, &perm);
2202:   *sectionNew = sNew;
2203:   return(0);
2204: }

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

2209:   Collective

2211:   Input Parameters:
2212: + sf - The SF
2213: - rootSection - Section defined on root space

2215:   Output Parameters:
2216: + remoteOffsets - root offsets in leaf storage, or NULL
2217: - leafSection - Section defined on the leaf space

2219:   Level: intermediate

2221: .seealso: PetscSFCreate()
2222: @*/
2223: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2224: {
2225:   PetscSF        embedSF;
2226:   const PetscInt *ilocal, *indices;
2227:   IS             selected;
2228:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

2232:   PetscSectionGetNumFields(rootSection, &numFields);
2233:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2234:   for (f = 0; f < numFields; ++f) {
2235:     PetscInt numComp = 0;
2236:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2237:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2238:   }
2239:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2240:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2241:   rpEnd = PetscMin(rpEnd,nroots);
2242:   rpEnd = PetscMax(rpStart,rpEnd);
2243:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2244:   ISGetIndices(selected, &indices);
2245:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2246:   ISRestoreIndices(selected, &indices);
2247:   ISDestroy(&selected);
2248:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2249:   if (nleaves && ilocal) {
2250:     for (i = 0; i < nleaves; ++i) {
2251:       lpStart = PetscMin(lpStart, ilocal[i]);
2252:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
2253:     }
2254:     ++lpEnd;
2255:   } else {
2256:     lpStart = 0;
2257:     lpEnd   = nleaves;
2258:   }
2259:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
2260:   /* Could fuse these at the cost of a copy and extra allocation */
2261:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2262:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2263:   for (f = 0; f < numFields; ++f) {
2264:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2265:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2266:   }
2267:   if (remoteOffsets) {
2268:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2269:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2270:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2271:   }
2272:   PetscSFDestroy(&embedSF);
2273:   PetscSectionSetUp(leafSection);
2274:   return(0);
2275: }

2277: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2278: {
2279:   PetscSF         embedSF;
2280:   const PetscInt *indices;
2281:   IS              selected;
2282:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2283:   PetscErrorCode  ierr;

2286:   *remoteOffsets = NULL;
2287:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2288:   if (numRoots < 0) return(0);
2289:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2290:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2291:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2292:   ISGetIndices(selected, &indices);
2293:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2294:   ISRestoreIndices(selected, &indices);
2295:   ISDestroy(&selected);
2296:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2297:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2298:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2299:   PetscSFDestroy(&embedSF);
2300:   return(0);
2301: }

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

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

2312:   Output Parameters:
2313: - sectionSF - The new SF

2315:   Note: Either rootSection or remoteOffsets can be specified

2317:   Level: intermediate

2319: .seealso: PetscSFCreate()
2320: @*/
2321: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2322: {
2323:   MPI_Comm          comm;
2324:   const PetscInt    *localPoints;
2325:   const PetscSFNode *remotePoints;
2326:   PetscInt          lpStart, lpEnd;
2327:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2328:   PetscInt          *localIndices;
2329:   PetscSFNode       *remoteIndices;
2330:   PetscInt          i, ind;
2331:   PetscErrorCode    ierr;

2339:   PetscObjectGetComm((PetscObject)sf,&comm);
2340:   PetscSFCreate(comm, sectionSF);
2341:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2342:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2343:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2344:   if (numRoots < 0) return(0);
2345:   for (i = 0; i < numPoints; ++i) {
2346:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2347:     PetscInt dof;

2349:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2350:       PetscSectionGetDof(leafSection, localPoint, &dof);
2351:       numIndices += dof;
2352:     }
2353:   }
2354:   PetscMalloc1(numIndices, &localIndices);
2355:   PetscMalloc1(numIndices, &remoteIndices);
2356:   /* Create new index graph */
2357:   for (i = 0, ind = 0; i < numPoints; ++i) {
2358:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2359:     PetscInt rank       = remotePoints[i].rank;

2361:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2362:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2363:       PetscInt loff, dof, d;

2365:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2366:       PetscSectionGetDof(leafSection, localPoint, &dof);
2367:       for (d = 0; d < dof; ++d, ++ind) {
2368:         localIndices[ind]        = loff+d;
2369:         remoteIndices[ind].rank  = rank;
2370:         remoteIndices[ind].index = remoteOffset+d;
2371:       }
2372:     }
2373:   }
2374:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2375:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2376:   return(0);
2377: }

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

2382:   Input Parameters:
2383: + section   - The PetscSection
2384: . obj       - A PetscObject which serves as the key for this index
2385: . clSection - Section giving the size of the closure of each point
2386: - clPoints  - IS giving the points in each closure

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

2390:   Level: intermediate

2392: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2393: @*/
2394: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2395: {

2399:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2400:   section->clObj     = obj;
2401:   PetscSectionDestroy(&section->clSection);
2402:   ISDestroy(&section->clPoints);
2403:   section->clSection = clSection;
2404:   section->clPoints  = clPoints;
2405:   return(0);
2406: }

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

2411:   Input Parameters:
2412: + section   - The PetscSection
2413: - obj       - A PetscObject which serves as the key for this index

2415:   Output Parameters:
2416: + clSection - Section giving the size of the closure of each point
2417: - clPoints  - IS giving the points in each closure

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

2421:   Level: intermediate

2423: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2424: @*/
2425: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2426: {
2428:   if (section->clObj == obj) {
2429:     if (clSection) *clSection = section->clSection;
2430:     if (clPoints)  *clPoints  = section->clPoints;
2431:   } else {
2432:     if (clSection) *clSection = NULL;
2433:     if (clPoints)  *clPoints  = NULL;
2434:   }
2435:   return(0);
2436: }

2438: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2439: {
2440:   PetscInt       i;

2444:   if (section->clObj != obj) {
2445:     PetscSectionDestroy(&section->clSection);
2446:     ISDestroy(&section->clPoints);
2447:   }
2448:   section->clObj  = obj;
2449:   PetscFree(section->clPerm);
2450:   PetscFree(section->clInvPerm);
2451:   section->clSize = clSize;
2452:   if (mode == PETSC_COPY_VALUES) {
2453:     PetscMalloc1(clSize, &section->clPerm);
2454:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2455:     PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2456:   } else if (mode == PETSC_OWN_POINTER) {
2457:     section->clPerm = clPerm;
2458:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2459:   PetscMalloc1(clSize, &section->clInvPerm);
2460:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2461:   return(0);
2462: }

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

2467:   Input Parameters:
2468: + section - The PetscSection
2469: . obj     - A PetscObject which serves as the key for this index
2470: - perm    - Permutation of the cell dof closure

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

2475:   Level: intermediate

2477: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2478: @*/
2479: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2480: {
2481:   const PetscInt *clPerm = NULL;
2482:   PetscInt        clSize = 0;
2483:   PetscErrorCode  ierr;

2486:   if (perm) {
2487:     ISGetLocalSize(perm, &clSize);
2488:     ISGetIndices(perm, &clPerm);
2489:   }
2490:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2491:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2492:   return(0);
2493: }

2495: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2496: {
2498:   if (section->clObj == obj) {
2499:     if (size) *size = section->clSize;
2500:     if (perm) *perm = section->clPerm;
2501:   } else {
2502:     if (size) *size = 0;
2503:     if (perm) *perm = NULL;
2504:   }
2505:   return(0);
2506: }

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

2511:   Input Parameters:
2512: + section   - The PetscSection
2513: - obj       - A PetscObject which serves as the key for this index

2515:   Output Parameter:
2516: . perm - The dof closure permutation

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

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

2523:   Level: intermediate

2525: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2526: @*/
2527: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2528: {
2529:   const PetscInt *clPerm;
2530:   PetscInt        clSize;
2531:   PetscErrorCode  ierr;

2534:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2535:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2536:   return(0);
2537: }

2539: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2540: {
2542:   if (section->clObj == obj) {
2543:     if (size) *size = section->clSize;
2544:     if (perm) *perm = section->clInvPerm;
2545:   } else {
2546:     if (size) *size = 0;
2547:     if (perm) *perm = NULL;
2548:   }
2549:   return(0);
2550: }

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

2555:   Input Parameters:
2556: + section   - The PetscSection
2557: - obj       - A PetscObject which serves as the key for this index

2559:   Output Parameters:
2560: + size - The dof closure size
2561: - perm - The dof closure permutation

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

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

2568:   Level: intermediate

2570: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2571: @*/
2572: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2573: {
2574:   const PetscInt *clPerm;
2575:   PetscInt        clSize;
2576:   PetscErrorCode  ierr;

2579:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2580:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2581:   return(0);
2582: }

2584: /*@
2585:   PetscSectionGetField - Get the subsection associated with a single field

2587:   Input Parameters:
2588: + s     - The PetscSection
2589: - field - The field number

2591:   Output Parameter:
2592: . subs  - The subsection for the given field

2594:   Level: intermediate

2596: .seealso: PetscSectionSetNumFields()
2597: @*/
2598: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2599: {
2603:   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);
2604:   *subs = s->field[field];
2605:   return(0);
2606: }

2608: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2609: PetscFunctionList PetscSectionSymList = NULL;

2611: /*@
2612:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2614:   Collective on MPI_Comm

2616:   Input Parameter:
2617: . comm - the MPI communicator

2619:   Output Parameter:
2620: . sym - pointer to the new set of symmetries

2622:   Level: developer

2624: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2625: @*/
2626: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2627: {

2632:   ISInitializePackage();
2633:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2634:   return(0);
2635: }

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

2640:   Collective on PetscSectionSym

2642:   Input Parameters:
2643: + sym    - The section symmetry object
2644: - method - The name of the section symmetry type

2646:   Level: developer

2648: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2649: @*/
2650: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2651: {
2652:   PetscErrorCode (*r)(PetscSectionSym);
2653:   PetscBool      match;

2658:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2659:   if (match) return(0);

2661:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2662:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2663:   if (sym->ops->destroy) {
2664:     (*sym->ops->destroy)(sym);
2665:     sym->ops->destroy = NULL;
2666:   }
2667:   (*r)(sym);
2668:   PetscObjectChangeTypeName((PetscObject)sym,method);
2669:   return(0);
2670: }


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

2676:   Not Collective

2678:   Input Parameter:
2679: . sym  - The section symmetry

2681:   Output Parameter:
2682: . type - The index set type name

2684:   Level: developer

2686: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2687: @*/
2688: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2689: {
2693:   *type = ((PetscObject)sym)->type_name;
2694:   return(0);
2695: }

2697: /*@C
2698:   PetscSectionSymRegister - Adds a new section symmetry implementation

2700:   Not Collective

2702:   Input Parameters:
2703: + name        - The name of a new user-defined creation routine
2704: - create_func - The creation routine itself

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

2709:   Level: developer

2711: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2712: @*/
2713: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2714: {

2718:   ISInitializePackage();
2719:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2720:   return(0);
2721: }

2723: /*@
2724:    PetscSectionSymDestroy - Destroys a section symmetry.

2726:    Collective on PetscSectionSym

2728:    Input Parameters:
2729: .  sym - the section symmetry

2731:    Level: developer

2733: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2734: @*/
2735: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2736: {
2737:   SymWorkLink    link,next;

2741:   if (!*sym) return(0);
2743:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2744:   if ((*sym)->ops->destroy) {
2745:     (*(*sym)->ops->destroy)(*sym);
2746:   }
2747:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2748:   for (link=(*sym)->workin; link; link=next) {
2749:     next = link->next;
2750:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2751:     PetscFree(link);
2752:   }
2753:   (*sym)->workin = NULL;
2754:   PetscHeaderDestroy(sym);
2755:   return(0);
2756: }

2758: /*@C
2759:    PetscSectionSymView - Displays a section symmetry

2761:    Collective on PetscSectionSym

2763:    Input Parameters:
2764: +  sym - the index set
2765: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2767:    Level: developer

2769: .seealso: PetscViewerASCIIOpen()
2770: @*/
2771: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2772: {

2777:   if (!viewer) {
2778:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2779:   }
2782:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2783:   if (sym->ops->view) {
2784:     (*sym->ops->view)(sym,viewer);
2785:   }
2786:   return(0);
2787: }

2789: /*@
2790:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2792:   Collective on PetscSection

2794:   Input Parameters:
2795: + section - the section describing data layout
2796: - sym - the symmetry describing the affect of orientation on the access of the data

2798:   Level: developer

2800: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2801: @*/
2802: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2803: {

2808:   PetscSectionSymDestroy(&(section->sym));
2809:   if (sym) {
2812:     PetscObjectReference((PetscObject) sym);
2813:   }
2814:   section->sym = sym;
2815:   return(0);
2816: }

2818: /*@
2819:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2821:   Not collective

2823:   Input Parameters:
2824: . section - the section describing data layout

2826:   Output Parameters:
2827: . sym - the symmetry describing the affect of orientation on the access of the data

2829:   Level: developer

2831: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2832: @*/
2833: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2834: {
2837:   *sym = section->sym;
2838:   return(0);
2839: }

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

2844:   Collective on PetscSection

2846:   Input Parameters:
2847: + section - the section describing data layout
2848: . field - the field number
2849: - sym - the symmetry describing the affect of orientation on the access of the data

2851:   Level: developer

2853: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2854: @*/
2855: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2856: {

2861:   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);
2862:   PetscSectionSetSym(section->field[field],sym);
2863:   return(0);
2864: }

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

2869:   Collective on PetscSection

2871:   Input Parameters:
2872: + section - the section describing data layout
2873: - field - the field number

2875:   Output Parameters:
2876: . sym - the symmetry describing the affect of orientation on the access of the data

2878:   Level: developer

2880: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2881: @*/
2882: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2883: {
2886:   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);
2887:   *sym = section->field[field]->sym;
2888:   return(0);
2889: }

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

2894:   Not collective

2896:   Input Parameters:
2897: + section - the section
2898: . numPoints - the number of points
2899: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2900:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2901:     context, see DMPlexGetConeOrientation()).

2903:   Output Parameter:
2904: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2905: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2906:     identity).

2908:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2909: .vb
2910:      const PetscInt    **perms;
2911:      const PetscScalar **rots;
2912:      PetscInt            lOffset;

2914:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2915:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2916:        PetscInt           point = points[2*i], dof, sOffset;
2917:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2918:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2920:        PetscSectionGetDof(section,point,&dof);
2921:        PetscSectionGetOffset(section,point,&sOffset);

2923:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2924:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2925:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2926:        lOffset += dof;
2927:      }
2928:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2929: .ve

2931:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2932: .vb
2933:      const PetscInt    **perms;
2934:      const PetscScalar **rots;
2935:      PetscInt            lOffset;

2937:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2938:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2939:        PetscInt           point = points[2*i], dof, sOffset;
2940:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2941:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2943:        PetscSectionGetDof(section,point,&dof);
2944:        PetscSectionGetOffset(section,point,&sOff);

2946:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2947:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2948:        offset += dof;
2949:      }
2950:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2951: .ve

2953:   Level: developer

2955: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2956: @*/
2957: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2958: {
2959:   PetscSectionSym sym;
2960:   PetscErrorCode  ierr;

2965:   if (perms) *perms = NULL;
2966:   if (rots)  *rots  = NULL;
2967:   sym = section->sym;
2968:   if (sym && (perms || rots)) {
2969:     SymWorkLink link;

2971:     if (sym->workin) {
2972:       link        = sym->workin;
2973:       sym->workin = sym->workin->next;
2974:     } else {
2975:       PetscNewLog(sym,&link);
2976:     }
2977:     if (numPoints > link->numPoints) {
2978:       PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2979:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2980:       link->numPoints = numPoints;
2981:     }
2982:     link->next   = sym->workout;
2983:     sym->workout = link;
2984:     PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2985:     PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2986:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2987:     if (perms) *perms = link->perms;
2988:     if (rots)  *rots  = link->rots;
2989:   }
2990:   return(0);
2991: }

2993: /*@C
2994:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

2996:   Not collective

2998:   Input Parameters:
2999: + section - the section
3000: . numPoints - the number of points
3001: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3002:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3003:     context, see DMPlexGetConeOrientation()).

3005:   Output Parameter:
3006: + perms - The permutations for the given orientations: set to NULL at conclusion
3007: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3009:   Level: developer

3011: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3012: @*/
3013: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3014: {
3015:   PetscSectionSym sym;

3019:   sym = section->sym;
3020:   if (sym && (perms || rots)) {
3021:     SymWorkLink *p,link;

3023:     for (p=&sym->workout; (link=*p); p=&link->next) {
3024:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3025:         *p          = link->next;
3026:         link->next  = sym->workin;
3027:         sym->workin = link;
3028:         if (perms) *perms = NULL;
3029:         if (rots)  *rots  = NULL;
3030:         return(0);
3031:       }
3032:     }
3033:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3034:   }
3035:   return(0);
3036: }

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

3041:   Not collective

3043:   Input Parameters:
3044: + section - the section
3045: . field - the field of the section
3046: . numPoints - the number of points
3047: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3048:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3049:     context, see DMPlexGetConeOrientation()).

3051:   Output Parameter:
3052: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3053: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3054:     identity).

3056:   Level: developer

3058: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3059: @*/
3060: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3061: {

3066:   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);
3067:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3068:   return(0);
3069: }

3071: /*@C
3072:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3074:   Not collective

3076:   Input Parameters:
3077: + section - the section
3078: . field - the field number
3079: . numPoints - the number of points
3080: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3081:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3082:     context, see DMPlexGetConeOrientation()).

3084:   Output Parameter:
3085: + perms - The permutations for the given orientations: set to NULL at conclusion
3086: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3088:   Level: developer

3090: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3091: @*/
3092: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3093: {

3098:   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);
3099:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3100:   return(0);
3101: }

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

3106:   Not collective

3108:   Input Parameter:
3109: . s - the global PetscSection

3111:   Output Parameters:
3112: . flg - the flag

3114:   Level: developer

3116: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3117: @*/
3118: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3119: {
3122:   *flg = s->useFieldOff;
3123:   return(0);
3124: }

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

3129:   Not collective

3131:   Input Parameters:
3132: + s   - the global PetscSection
3133: - flg - the flag

3135:   Level: developer

3137: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3138: @*/
3139: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3140: {
3143:   s->useFieldOff = flg;
3144:   return(0);
3145: }

3147: #define PetscSectionExpandPoints_Loop(TYPE) \
3148: { \
3149:   PetscInt i, n, o0, o1, size; \
3150:   TYPE *a0 = (TYPE*)origArray, *a1; \
3151:   PetscSectionGetStorageSize(s, &size); \
3152:   PetscMalloc1(size, &a1); \
3153:   for (i=0; i<npoints; i++) { \
3154:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3155:     PetscSectionGetOffset(s, i, &o1); \
3156:     PetscSectionGetDof(s, i, &n); \
3157:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3158:   } \
3159:   *newArray = (void*)a1; \
3160: }

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

3165:   Not collective

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

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

3177:   Level: developer

3179: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3180: @*/
3181: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3182: {
3183:   PetscSection        s;
3184:   const PetscInt      *points_;
3185:   PetscInt            i, n, npoints, pStart, pEnd;
3186:   PetscMPIInt         unitsize;
3187:   PetscErrorCode      ierr;

3195:   MPI_Type_size(dataType, &unitsize);
3196:   ISGetLocalSize(points, &npoints);
3197:   ISGetIndices(points, &points_);
3198:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3199:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3200:   PetscSectionSetChart(s, 0, npoints);
3201:   for (i=0; i<npoints; i++) {
3202:     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);
3203:     PetscSectionGetDof(origSection, points_[i], &n);
3204:     PetscSectionSetDof(s, i, n);
3205:   }
3206:   PetscSectionSetUp(s);
3207:   if (newArray) {
3208:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3209:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3210:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3211:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3212:   }
3213:   if (newSection) {
3214:     *newSection = s;
3215:   } else {
3216:     PetscSectionDestroy(&s);
3217:   }
3218:   return(0);
3219: }