Actual source code: vsectionis.c

petsc-master 2019-08-19
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

 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

 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

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:     PetscArraycmp(idx1, idx2, ncdof, 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:       PetscArraycmp(idx1, idx2, nfcdof, 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:   PetscArrayzero(s->atlasDof, pEnd - pStart);
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 defined(PETSC_USE_DEBUG)
727:   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);
728: #endif
729:   s->atlasDof[point - s->pStart] += numDof;
730:   return(0);
731: }

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

736:   Not collective

738:   Input Parameters:
739: + s - the PetscSection
740: . point - the point
741: - field - the field

743:   Output Parameter:
744: . numDof - the number of dof

746:   Level: intermediate

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

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

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

764:   Not collective

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

772:   Level: intermediate

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

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

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

790:   Not collective

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

798:   Level: intermediate

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

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

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

816:   Not collective

818:   Input Parameters:
819: + s - the PetscSection
820: - point - the point

822:   Output Parameter:
823: . numDof - the number of dof which are fixed by constraints

825:   Level: intermediate

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

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

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

845:   Not collective

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

852:   Level: intermediate

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

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

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

872:   Not collective

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

879:   Level: intermediate

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

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

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

899:   Not collective

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

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

909:   Level: intermediate

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

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

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

928:   Not collective

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

936:   Level: intermediate

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

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

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

954:   Not collective

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

962:   Level: intermediate

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

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

977: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
978: {

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

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

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

995:   Not collective

997:   Input Parameter:
998: . s - the PetscSection

1000:   Level: intermediate

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

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

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

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

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

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

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

1061:   Not collective

1063:   Input Parameters:
1064: . s - the PetscSection

1066:   Output Parameter:
1067: . maxDof - the maximum dof

1069:   Level: intermediate

1071: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1072: @*/
1073: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1074: {
1078:   *maxDof = s->maxDof;
1079:   return(0);
1080: }

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

1085:   Not collective

1087:   Input Parameter:
1088: . s - the PetscSection

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

1093:   Level: intermediate

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

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

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

1112:   Not collective

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

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

1121:   Level: intermediate

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

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

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

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

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

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

1155:   Level: developer

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

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

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

1241: /*@
1242:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1243:   the local section and an SF describing the section point overlap.

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

1252:   Output Parameter:
1253:   . gsection - The PetscSection for the global field layout

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

1257:   Level: developer

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

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

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

1338: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1339: {
1340:   PetscInt       pStart, pEnd, p, localSize = 0;

1344:   PetscSectionGetChart(s, &pStart, &pEnd);
1345:   for (p = pStart; p < pEnd; ++p) {
1346:     PetscInt dof;

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

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

1361:   Input Parameters:
1362: + comm - The MPI_Comm
1363: - s    - The PetscSection

1365:   Output Parameter:
1366: . layout - The layout for the section

1368:   Level: developer

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

1380:   PetscSectionGetChart(s, &pStart, &pEnd);
1381:   for (p = pStart; p < pEnd; ++p) {
1382:     PetscInt dof,cdof;

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

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

1398:   Not collective

1400:   Input Parameters:
1401: + s - the PetscSection
1402: - point - the point

1404:   Output Parameter:
1405: . offset - the offset

1407:   Level: intermediate

1409: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1410: @*/
1411: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1412: {
1416: #if defined(PETSC_USE_DEBUG)
1417:   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);
1418: #endif
1419:   *offset = s->atlasOff[point - s->pStart];
1420:   return(0);
1421: }

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

1426:   Not collective

1428:   Input Parameters:
1429: + s - the PetscSection
1430: . point - the point
1431: - offset - the offset

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

1435:   Level: intermediate

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

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

1451:   Not collective

1453:   Input Parameters:
1454: + s - the PetscSection
1455: . point - the point
1456: - field - the field

1458:   Output Parameter:
1459: . offset - the offset

1461:   Level: intermediate

1463: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1464: @*/
1465: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1466: {

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

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

1480:   Not collective

1482:   Input Parameters:
1483: + s - the PetscSection
1484: . point - the point
1485: . field - the field
1486: - offset - the offset

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

1490:   Level: intermediate

1492: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1493: @*/
1494: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1495: {

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

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

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

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

1524:   Not collective

1526:   Input Parameter:
1527: . s - the PetscSection

1529:   Output Parameters:
1530: + start - the minimum offset
1531: - end   - one more than the maximum offset

1533:   Level: intermediate

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

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

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

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

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

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

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

1603:     PetscMalloc1(maxCdof, &indices);
1604:     for (p = pStart; p < pEnd; ++p) {
1605:       PetscInt cdof;

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

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

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

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

1643:   if (!len) return(0);
1644:   for (i = 0; i < len; ++i) {
1645:     PetscInt nf, pStarti, pEndi;

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

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

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

1673:     for (i = 0, f = 0; i < len; ++i) {
1674:       PetscInt nf, fi, pStarti, pEndi;
1675:       PetscInt fdof = 0, cfdof = 0;

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

1697:     PetscMalloc1(maxCdof, &indices);
1698:     for (p = pStart; p < pEnd; ++p) {
1699:       PetscInt cdof;

1701:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1702:       if (cdof) {
1703:         PetscInt dof, numConst = 0, fOff = 0;

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

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

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

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

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

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

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

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

1808: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1809: {
1810:   PetscInt       p;
1811:   PetscMPIInt    rank;

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

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

1841: /*@C
1842:   PetscSectionView - Views a PetscSection

1844:   Collective on PetscSection

1846:   Input Parameters:
1847: + s - the PetscSection object to view
1848: - v - the viewer

1850:   Level: developer

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

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

1880: /*@
1881:   PetscSectionReset - Frees all section data.

1883:   Not collective

1885:   Input Parameters:
1886: . s - the PetscSection

1888:   Level: developer

1890: .seealso: PetscSection, PetscSectionCreate()
1891: @*/
1892: PetscErrorCode PetscSectionReset(PetscSection s)
1893: {
1894:   PetscInt       f;

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

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

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

1928:   Not collective

1930:   Input Parameters:
1931: . s - the PetscSection

1933:   Level: developer

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

1938: .seealso: PetscSection, PetscSectionCreate()
1939: @*/
1940: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1941: {

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

1956: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1957: {
1958:   const PetscInt p = point - s->pStart;

1962:   *values = &baseArray[s->atlasOff[p]];
1963:   return(0);
1964: }

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

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

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

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

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

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

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

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

2042: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2043: {
2047:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2048:   return(0);
2049: }

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

2054:   Input Parameters:
2055: + s     - The PetscSection
2056: - point - The point

2058:   Output Parameter:
2059: . indices - The constrained dofs

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

2063:   Level: advanced

2065: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2066: @*/
2067: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2068: {

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

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

2082:   Input Parameters:
2083: + s     - The PetscSection
2084: . point - The point
2085: - indices - The constrained dofs

2087:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2089:   Level: advanced

2091: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2092: @*/
2093: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2094: {

2099:   if (s->bc) {
2100:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2101:   }
2102:   return(0);
2103: }

2105: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2106: {

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

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

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

2127: /*@
2128:   PetscSectionPermute - Reorder the section according to the input point permutation

2130:   Collective on PetscSection

2132:   Input Parameter:
2133: + section - The PetscSection object
2134: - perm - The point permutation, old point p becomes new point perm[p]

2136:   Output Parameter:
2137: . sectionNew - The permuted PetscSection

2139:   Level: intermediate

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

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

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

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

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

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

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

2214:   Collective

2216:   Input Parameters:
2217: + sf - The SF
2218: - rootSection - Section defined on root space

2220:   Output Parameters:
2221: + remoteOffsets - root offsets in leaf storage, or NULL
2222: - leafSection - Section defined on the leaf space

2224:   Level: intermediate

2226: .seealso: PetscSFCreate()
2227: @*/
2228: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2229: {
2230:   PetscSF        embedSF;
2231:   const PetscInt *indices;
2232:   IS             selected;
2233:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f;
2234:   PetscBool      sub,lsub;

2238:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2239:   PetscSectionGetNumFields(rootSection, &numFields);
2240:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2241:   for (f = 0; f < numFields; ++f) {
2242:     PetscInt numComp = 0;
2243:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2244:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2245:   }
2246:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2247:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2248:   rpEnd = PetscMin(rpEnd,nroots);
2249:   rpEnd = PetscMax(rpStart,rpEnd);
2250:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2251:   lsub = (PetscBool)(nroots != rpEnd - rpStart);
2252:   MPIU_Allreduce(&lsub,&sub,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)sf));
2253:   if (sub) {
2254:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2255:     ISGetIndices(selected, &indices);
2256:     PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2257:     ISRestoreIndices(selected, &indices);
2258:     ISDestroy(&selected);
2259:   } else {
2260:     PetscObjectReference((PetscObject)sf);
2261:     embedSF = sf;
2262:   }
2263:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2264:   lpEnd++;

2266:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
2267:   /* Could fuse these at the cost of a copy and extra allocation */
2268:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2269:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2270:   for (f = 0; f < numFields; ++f) {
2271:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2272:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2273:   }
2274:   if (remoteOffsets) {
2275:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2276:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2277:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2278:   }
2279:   PetscSFDestroy(&embedSF);
2280:   PetscSectionSetUp(leafSection);
2281:   PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2282:   return(0);
2283: }

2285: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2286: {
2287:   PetscSF         embedSF;
2288:   const PetscInt *indices;
2289:   IS              selected;
2290:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2291:   PetscErrorCode  ierr;

2294:   *remoteOffsets = NULL;
2295:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2296:   if (numRoots < 0) return(0);
2297:   PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2298:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2299:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2300:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2301:   ISGetIndices(selected, &indices);
2302:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2303:   ISRestoreIndices(selected, &indices);
2304:   ISDestroy(&selected);
2305:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2306:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2307:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2308:   PetscSFDestroy(&embedSF);
2309:   PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2310:   return(0);
2311: }

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

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

2322:   Output Parameters:
2323: - sectionSF - The new SF

2325:   Note: Either rootSection or remoteOffsets can be specified

2327:   Level: intermediate

2329: .seealso: PetscSFCreate()
2330: @*/
2331: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2332: {
2333:   MPI_Comm          comm;
2334:   const PetscInt    *localPoints;
2335:   const PetscSFNode *remotePoints;
2336:   PetscInt          lpStart, lpEnd;
2337:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2338:   PetscInt          *localIndices;
2339:   PetscSFNode       *remoteIndices;
2340:   PetscInt          i, ind;
2341:   PetscErrorCode    ierr;

2349:   PetscObjectGetComm((PetscObject)sf,&comm);
2350:   PetscSFCreate(comm, sectionSF);
2351:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2352:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2353:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2354:   if (numRoots < 0) return(0);
2355:   PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2356:   for (i = 0; i < numPoints; ++i) {
2357:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2358:     PetscInt dof;

2360:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2361:       PetscSectionGetDof(leafSection, localPoint, &dof);
2362:       numIndices += dof;
2363:     }
2364:   }
2365:   PetscMalloc1(numIndices, &localIndices);
2366:   PetscMalloc1(numIndices, &remoteIndices);
2367:   /* Create new index graph */
2368:   for (i = 0, ind = 0; i < numPoints; ++i) {
2369:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2370:     PetscInt rank       = remotePoints[i].rank;

2372:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2373:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2374:       PetscInt loff, dof, d;

2376:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2377:       PetscSectionGetDof(leafSection, localPoint, &dof);
2378:       for (d = 0; d < dof; ++d, ++ind) {
2379:         localIndices[ind]        = loff+d;
2380:         remoteIndices[ind].rank  = rank;
2381:         remoteIndices[ind].index = remoteOffset+d;
2382:       }
2383:     }
2384:   }
2385:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2386:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2387:   PetscSFSetUp(*sectionSF);
2388:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2389:   return(0);
2390: }

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

2395:   Input Parameters:
2396: + section   - The PetscSection
2397: . obj       - A PetscObject which serves as the key for this index
2398: . clSection - Section giving the size of the closure of each point
2399: - clPoints  - IS giving the points in each closure

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

2403:   Level: intermediate

2405: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2406: @*/
2407: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2408: {

2412:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2413:   section->clObj     = obj;
2414:   PetscSectionDestroy(&section->clSection);
2415:   ISDestroy(&section->clPoints);
2416:   section->clSection = clSection;
2417:   section->clPoints  = clPoints;
2418:   return(0);
2419: }

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

2424:   Input Parameters:
2425: + section   - The PetscSection
2426: - obj       - A PetscObject which serves as the key for this index

2428:   Output Parameters:
2429: + clSection - Section giving the size of the closure of each point
2430: - clPoints  - IS giving the points in each closure

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

2434:   Level: intermediate

2436: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2437: @*/
2438: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2439: {
2441:   if (section->clObj == obj) {
2442:     if (clSection) *clSection = section->clSection;
2443:     if (clPoints)  *clPoints  = section->clPoints;
2444:   } else {
2445:     if (clSection) *clSection = NULL;
2446:     if (clPoints)  *clPoints  = NULL;
2447:   }
2448:   return(0);
2449: }

2451: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2452: {
2453:   PetscInt       i;

2457:   if (section->clObj != obj) {
2458:     PetscSectionDestroy(&section->clSection);
2459:     ISDestroy(&section->clPoints);
2460:   }
2461:   section->clObj  = obj;
2462:   PetscFree(section->clPerm);
2463:   PetscFree(section->clInvPerm);
2464:   section->clSize = clSize;
2465:   if (mode == PETSC_COPY_VALUES) {
2466:     PetscMalloc1(clSize, &section->clPerm);
2467:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2468:     PetscArraycpy(section->clPerm, clPerm, clSize);
2469:   } else if (mode == PETSC_OWN_POINTER) {
2470:     section->clPerm = clPerm;
2471:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2472:   PetscMalloc1(clSize, &section->clInvPerm);
2473:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2474:   return(0);
2475: }

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

2480:   Input Parameters:
2481: + section - The PetscSection
2482: . obj     - A PetscObject which serves as the key for this index
2483: - perm    - Permutation of the cell dof closure

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

2488:   Level: intermediate

2490: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2491: @*/
2492: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2493: {
2494:   const PetscInt *clPerm = NULL;
2495:   PetscInt        clSize = 0;
2496:   PetscErrorCode  ierr;

2499:   if (perm) {
2500:     ISGetLocalSize(perm, &clSize);
2501:     ISGetIndices(perm, &clPerm);
2502:   }
2503:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2504:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2505:   return(0);
2506: }

2508: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2509: {
2511:   if (section->clObj == obj) {
2512:     if (size) *size = section->clSize;
2513:     if (perm) *perm = section->clPerm;
2514:   } else {
2515:     if (size) *size = 0;
2516:     if (perm) *perm = NULL;
2517:   }
2518:   return(0);
2519: }

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

2524:   Input Parameters:
2525: + section   - The PetscSection
2526: - obj       - A PetscObject which serves as the key for this index

2528:   Output Parameter:
2529: . perm - The dof closure permutation

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

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

2536:   Level: intermediate

2538: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2539: @*/
2540: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2541: {
2542:   const PetscInt *clPerm;
2543:   PetscInt        clSize;
2544:   PetscErrorCode  ierr;

2547:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2548:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2549:   return(0);
2550: }

2552: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2553: {
2555:   if (section->clObj == obj) {
2556:     if (size) *size = section->clSize;
2557:     if (perm) *perm = section->clInvPerm;
2558:   } else {
2559:     if (size) *size = 0;
2560:     if (perm) *perm = NULL;
2561:   }
2562:   return(0);
2563: }

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

2568:   Input Parameters:
2569: + section   - The PetscSection
2570: - obj       - A PetscObject which serves as the key for this index

2572:   Output Parameters:
2573: + size - The dof closure size
2574: - perm - The dof closure permutation

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

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

2581:   Level: intermediate

2583: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2584: @*/
2585: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2586: {
2587:   const PetscInt *clPerm;
2588:   PetscInt        clSize;
2589:   PetscErrorCode  ierr;

2592:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2593:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2594:   return(0);
2595: }

2597: /*@
2598:   PetscSectionGetField - Get the subsection associated with a single field

2600:   Input Parameters:
2601: + s     - The PetscSection
2602: - field - The field number

2604:   Output Parameter:
2605: . subs  - The subsection for the given field

2607:   Level: intermediate

2609: .seealso: PetscSectionSetNumFields()
2610: @*/
2611: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2612: {
2616:   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);
2617:   *subs = s->field[field];
2618:   return(0);
2619: }

2621: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2622: PetscFunctionList PetscSectionSymList = NULL;

2624: /*@
2625:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2627:   Collective

2629:   Input Parameter:
2630: . comm - the MPI communicator

2632:   Output Parameter:
2633: . sym - pointer to the new set of symmetries

2635:   Level: developer

2637: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2638: @*/
2639: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2640: {

2645:   ISInitializePackage();
2646:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2647:   return(0);
2648: }

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

2653:   Collective on PetscSectionSym

2655:   Input Parameters:
2656: + sym    - The section symmetry object
2657: - method - The name of the section symmetry type

2659:   Level: developer

2661: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2662: @*/
2663: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2664: {
2665:   PetscErrorCode (*r)(PetscSectionSym);
2666:   PetscBool      match;

2671:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2672:   if (match) return(0);

2674:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2675:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2676:   if (sym->ops->destroy) {
2677:     (*sym->ops->destroy)(sym);
2678:     sym->ops->destroy = NULL;
2679:   }
2680:   (*r)(sym);
2681:   PetscObjectChangeTypeName((PetscObject)sym,method);
2682:   return(0);
2683: }


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

2689:   Not Collective

2691:   Input Parameter:
2692: . sym  - The section symmetry

2694:   Output Parameter:
2695: . type - The index set type name

2697:   Level: developer

2699: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2700: @*/
2701: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2702: {
2706:   *type = ((PetscObject)sym)->type_name;
2707:   return(0);
2708: }

2710: /*@C
2711:   PetscSectionSymRegister - Adds a new section symmetry implementation

2713:   Not Collective

2715:   Input Parameters:
2716: + name        - The name of a new user-defined creation routine
2717: - create_func - The creation routine itself

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

2722:   Level: developer

2724: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2725: @*/
2726: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2727: {

2731:   ISInitializePackage();
2732:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2733:   return(0);
2734: }

2736: /*@
2737:    PetscSectionSymDestroy - Destroys a section symmetry.

2739:    Collective on PetscSectionSym

2741:    Input Parameters:
2742: .  sym - the section symmetry

2744:    Level: developer

2746: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2747: @*/
2748: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2749: {
2750:   SymWorkLink    link,next;

2754:   if (!*sym) return(0);
2756:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2757:   if ((*sym)->ops->destroy) {
2758:     (*(*sym)->ops->destroy)(*sym);
2759:   }
2760:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2761:   for (link=(*sym)->workin; link; link=next) {
2762:     next = link->next;
2763:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2764:     PetscFree(link);
2765:   }
2766:   (*sym)->workin = NULL;
2767:   PetscHeaderDestroy(sym);
2768:   return(0);
2769: }

2771: /*@C
2772:    PetscSectionSymView - Displays a section symmetry

2774:    Collective on PetscSectionSym

2776:    Input Parameters:
2777: +  sym - the index set
2778: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2780:    Level: developer

2782: .seealso: PetscViewerASCIIOpen()
2783: @*/
2784: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2785: {

2790:   if (!viewer) {
2791:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2792:   }
2795:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2796:   if (sym->ops->view) {
2797:     (*sym->ops->view)(sym,viewer);
2798:   }
2799:   return(0);
2800: }

2802: /*@
2803:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2805:   Collective on PetscSection

2807:   Input Parameters:
2808: + section - the section describing data layout
2809: - sym - the symmetry describing the affect of orientation on the access of the data

2811:   Level: developer

2813: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2814: @*/
2815: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2816: {

2821:   PetscSectionSymDestroy(&(section->sym));
2822:   if (sym) {
2825:     PetscObjectReference((PetscObject) sym);
2826:   }
2827:   section->sym = sym;
2828:   return(0);
2829: }

2831: /*@
2832:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2834:   Not collective

2836:   Input Parameters:
2837: . section - the section describing data layout

2839:   Output Parameters:
2840: . sym - the symmetry describing the affect of orientation on the access of the data

2842:   Level: developer

2844: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2845: @*/
2846: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2847: {
2850:   *sym = section->sym;
2851:   return(0);
2852: }

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

2857:   Collective on PetscSection

2859:   Input Parameters:
2860: + section - the section describing data layout
2861: . field - the field number
2862: - sym - the symmetry describing the affect of orientation on the access of the data

2864:   Level: developer

2866: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2867: @*/
2868: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2869: {

2874:   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);
2875:   PetscSectionSetSym(section->field[field],sym);
2876:   return(0);
2877: }

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

2882:   Collective on PetscSection

2884:   Input Parameters:
2885: + section - the section describing data layout
2886: - field - the field number

2888:   Output Parameters:
2889: . sym - the symmetry describing the affect of orientation on the access of the data

2891:   Level: developer

2893: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2894: @*/
2895: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2896: {
2899:   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);
2900:   *sym = section->field[field]->sym;
2901:   return(0);
2902: }

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

2907:   Not collective

2909:   Input Parameters:
2910: + section - the section
2911: . numPoints - the number of points
2912: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2913:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2914:     context, see DMPlexGetConeOrientation()).

2916:   Output Parameter:
2917: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2918: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2919:     identity).

2921:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2922: .vb
2923:      const PetscInt    **perms;
2924:      const PetscScalar **rots;
2925:      PetscInt            lOffset;

2927:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2928:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2929:        PetscInt           point = points[2*i], dof, sOffset;
2930:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2931:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2933:        PetscSectionGetDof(section,point,&dof);
2934:        PetscSectionGetOffset(section,point,&sOffset);

2936:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2937:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2938:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2939:        lOffset += dof;
2940:      }
2941:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2942: .ve

2944:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2945: .vb
2946:      const PetscInt    **perms;
2947:      const PetscScalar **rots;
2948:      PetscInt            lOffset;

2950:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2951:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2952:        PetscInt           point = points[2*i], dof, sOffset;
2953:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2954:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2956:        PetscSectionGetDof(section,point,&dof);
2957:        PetscSectionGetOffset(section,point,&sOff);

2959:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2960:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2961:        offset += dof;
2962:      }
2963:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2964: .ve

2966:   Level: developer

2968: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2969: @*/
2970: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2971: {
2972:   PetscSectionSym sym;
2973:   PetscErrorCode  ierr;

2978:   if (perms) *perms = NULL;
2979:   if (rots)  *rots  = NULL;
2980:   sym = section->sym;
2981:   if (sym && (perms || rots)) {
2982:     SymWorkLink link;

2984:     if (sym->workin) {
2985:       link        = sym->workin;
2986:       sym->workin = sym->workin->next;
2987:     } else {
2988:       PetscNewLog(sym,&link);
2989:     }
2990:     if (numPoints > link->numPoints) {
2991:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
2992:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2993:       link->numPoints = numPoints;
2994:     }
2995:     link->next   = sym->workout;
2996:     sym->workout = link;
2997:     PetscArrayzero((PetscInt**)link->perms,numPoints);
2998:     PetscArrayzero((PetscInt**)link->rots,numPoints);
2999:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3000:     if (perms) *perms = link->perms;
3001:     if (rots)  *rots  = link->rots;
3002:   }
3003:   return(0);
3004: }

3006: /*@C
3007:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3009:   Not collective

3011:   Input Parameters:
3012: + section - the section
3013: . numPoints - the number of points
3014: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3015:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3016:     context, see DMPlexGetConeOrientation()).

3018:   Output Parameter:
3019: + perms - The permutations for the given orientations: set to NULL at conclusion
3020: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3022:   Level: developer

3024: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3025: @*/
3026: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3027: {
3028:   PetscSectionSym sym;

3032:   sym = section->sym;
3033:   if (sym && (perms || rots)) {
3034:     SymWorkLink *p,link;

3036:     for (p=&sym->workout; (link=*p); p=&link->next) {
3037:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3038:         *p          = link->next;
3039:         link->next  = sym->workin;
3040:         sym->workin = link;
3041:         if (perms) *perms = NULL;
3042:         if (rots)  *rots  = NULL;
3043:         return(0);
3044:       }
3045:     }
3046:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3047:   }
3048:   return(0);
3049: }

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

3054:   Not collective

3056:   Input Parameters:
3057: + section - the section
3058: . field - the field of the section
3059: . numPoints - the number of points
3060: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3061:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3062:     context, see DMPlexGetConeOrientation()).

3064:   Output Parameter:
3065: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3066: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3067:     identity).

3069:   Level: developer

3071: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3072: @*/
3073: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3074: {

3079:   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);
3080:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3081:   return(0);
3082: }

3084: /*@C
3085:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3087:   Not collective

3089:   Input Parameters:
3090: + section - the section
3091: . field - the field number
3092: . numPoints - the number of points
3093: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3094:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3095:     context, see DMPlexGetConeOrientation()).

3097:   Output Parameter:
3098: + perms - The permutations for the given orientations: set to NULL at conclusion
3099: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3101:   Level: developer

3103: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3104: @*/
3105: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3106: {

3111:   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);
3112:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3113:   return(0);
3114: }

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

3119:   Not collective

3121:   Input Parameter:
3122: . s - the global PetscSection

3124:   Output Parameters:
3125: . flg - the flag

3127:   Level: developer

3129: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3130: @*/
3131: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3132: {
3135:   *flg = s->useFieldOff;
3136:   return(0);
3137: }

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

3142:   Not collective

3144:   Input Parameters:
3145: + s   - the global PetscSection
3146: - flg - the flag

3148:   Level: developer

3150: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3151: @*/
3152: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3153: {
3156:   s->useFieldOff = flg;
3157:   return(0);
3158: }

3160: #define PetscSectionExpandPoints_Loop(TYPE) \
3161: { \
3162:   PetscInt i, n, o0, o1, size; \
3163:   TYPE *a0 = (TYPE*)origArray, *a1; \
3164:   PetscSectionGetStorageSize(s, &size); \
3165:   PetscMalloc1(size, &a1); \
3166:   for (i=0; i<npoints; i++) { \
3167:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3168:     PetscSectionGetOffset(s, i, &o1); \
3169:     PetscSectionGetDof(s, i, &n); \
3170:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3171:   } \
3172:   *newArray = (void*)a1; \
3173: }

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

3178:   Not collective

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

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

3190:   Level: developer

3192: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3193: @*/
3194: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3195: {
3196:   PetscSection        s;
3197:   const PetscInt      *points_;
3198:   PetscInt            i, n, npoints, pStart, pEnd;
3199:   PetscMPIInt         unitsize;
3200:   PetscErrorCode      ierr;

3208:   MPI_Type_size(dataType, &unitsize);
3209:   ISGetLocalSize(points, &npoints);
3210:   ISGetIndices(points, &points_);
3211:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3212:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3213:   PetscSectionSetChart(s, 0, npoints);
3214:   for (i=0; i<npoints; i++) {
3215:     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);
3216:     PetscSectionGetDof(origSection, points_[i], &n);
3217:     PetscSectionSetDof(s, i, n);
3218:   }
3219:   PetscSectionSetUp(s);
3220:   if (newArray) {
3221:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3222:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3223:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3224:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3225:   }
3226:   if (newSection) {
3227:     *newSection = s;
3228:   } else {
3229:     PetscSectionDestroy(&s);
3230:   }
3231:   return(0);
3232: }