Actual source code: section.c

petsc-master 2019-12-03
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

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

 10: PetscClassId PETSC_SECTION_CLASSID;

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

 15:   Collective

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

 21:   Level: beginner

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

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

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

 44:   ISInitializePackage();

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

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

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

 74:   Collective

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

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

 82:   Level: intermediate

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

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

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

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

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

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

158:   Collective

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

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

166:   Level: beginner

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

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

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

185:   Collective on PetscSection

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

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

193:   Level: intermediate

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

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

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

215:   Collective on PetscSection

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

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

224:   Level: intermediate

226:   Notes:
227:   Field names are disregarded.

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

244:   flg = PETSC_FALSE;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

320:   Not collective

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

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

328:   Level: intermediate

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

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

344:   Not collective

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

350:   Level: intermediate

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

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

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

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

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

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

383:   Not Collective

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

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

392:   Level: intermediate

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

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

409:   Not Collective

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

416:   Level: intermediate

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

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

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

436:   Not collective

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

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

445:   Level: intermediate

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

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

462:   Not collective

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

469:   Level: intermediate

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

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

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

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

497:   Not collective

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

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

506:   Level: intermediate

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

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

522:   Not collective

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

529:   Level: intermediate

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

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

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

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

559:   Not collective

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

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

567:   Level: intermediate

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

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

582:   Not collective

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

588:   Level: intermediate

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

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

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

613:   Not collective

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

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

621:   Level: intermediate

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

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

637:   Not collective

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

643:   Not collective

645:   Level: intermediate

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

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

661:   Not collective

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

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

670:   Level: intermediate

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

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

689:   Not collective

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

696:   Level: intermediate

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

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

712:   Not collective

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

719:   Level: intermediate

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

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

737:   Not collective

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

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

747:   Level: intermediate

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

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

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

765:   Not collective

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

773:   Level: intermediate

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

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

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

791:   Not collective

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

799:   Level: intermediate

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

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

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

817:   Not collective

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

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

826:   Level: intermediate

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

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

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

846:   Not collective

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

853:   Level: intermediate

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

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

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

873:   Not collective

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

880:   Level: intermediate

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

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

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

900:   Not collective

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

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

910:   Level: intermediate

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

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

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

929:   Not collective

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

937:   Level: intermediate

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

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

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

955:   Not collective

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

963:   Level: intermediate

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

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

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

981:   Not collective

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

986:   Level: advanced

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

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

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

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

1008:   Not collective

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

1013:   Level: intermediate

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

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

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

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

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

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

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

1074:   Not collective

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

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

1082:   Level: intermediate

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

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

1098:   Not collective

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

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

1106:   Level: intermediate

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

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

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

1125:   Not collective

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

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

1134:   Level: intermediate

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

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

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

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

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

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

1168:   Level: intermediate

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

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

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

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

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

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

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

1270:   Level: advanced

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

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

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

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

1354:   Collective on comm

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

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

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

1365:   Level: advanced

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

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

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

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

1392:   Collective on comm

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

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

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

1403:   Level: advanced

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

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

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

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

1433:   Not collective

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

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

1442:   Level: intermediate

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

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

1461:   Not collective

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

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

1470:   Level: intermediate

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

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

1486:   Not collective

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

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

1496:   Level: intermediate

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

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

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

1515:   Not collective

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

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

1525:   Level: intermediate

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

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

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

1543:   Not collective

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

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

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

1556:   Level: advanced

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

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

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

1578:   Not collective

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

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

1587:   Level: intermediate

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

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

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

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

1616:   Collective on s

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

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

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

1628:   Level: advanced

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

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

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

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

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

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

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

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

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

1713:   Collective on s

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

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

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

1724:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

1825:   Collective on s

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

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

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

1836:   Level: advanced

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

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

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

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

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

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

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

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

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

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

1953:    Collective on PetscSection

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

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

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

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

1976:   Collective on PetscSection

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

1982:   Level: beginner

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

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

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

2015:   Not collective

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

2020:   Level: beginner

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

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

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

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

2060:   Not collective

2062:   Input Parameters:
2063: . s - the PetscSection

2065:   Level: beginner

2067: .seealso: PetscSection, PetscSectionCreate()
2068: @*/
2069: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2070: {

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

2085: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2086: {
2087:   const PetscInt p = point - s->pStart;

2091:   *values = &baseArray[s->atlasOff[p]];
2092:   return(0);
2093: }

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

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

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

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

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

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

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

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

2171: /*@C
2172:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2174:   Not collective

2176:   Input Parameter:
2177: . s - The PetscSection

2179:   Output Parameter:
2180: . hasConstraints - flag indicating that the section has constrained dofs

2182:   Level: intermediate

2184: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2185: @*/
2186: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2187: {
2191:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2192:   return(0);
2193: }

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

2198:   Not collective

2200:   Input Parameters:
2201: + s     - The PetscSection
2202: - point - The point

2204:   Output Parameter:
2205: . indices - The constrained dofs

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

2209:   Level: intermediate

2211: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2212: @*/
2213: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2214: {

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

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

2228:   Not collective

2230:   Input Parameters:
2231: + s     - The PetscSection
2232: . point - The point
2233: - indices - The constrained dofs

2235:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2237:   Level: intermediate

2239: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2240: @*/
2241: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2242: {

2247:   if (s->bc) {
2248:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2249:   }
2250:   return(0);
2251: }

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

2256:   Not collective

2258:   Input Parameters:
2259: + s     - The PetscSection
2260: . field  - The field number
2261: - point - The point

2263:   Output Parameter:
2264: . indices - The constrained dofs

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

2268:   Level: intermediate

2270: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2271: @*/
2272: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2273: {

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

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

2286:   Not collective

2288:   Input Parameters:
2289: + s       - The PetscSection
2290: . point   - The point
2291: . field   - The field number
2292: - indices - The constrained dofs

2294:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2296:   Level: intermediate

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

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

2311: /*@
2312:   PetscSectionPermute - Reorder the section according to the input point permutation

2314:   Collective on PetscSection

2316:   Input Parameter:
2317: + section - The PetscSection object
2318: - perm - The point permutation, old point p becomes new point perm[p]

2320:   Output Parameter:
2321: . sectionNew - The permuted PetscSection

2323:   Level: intermediate

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

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

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

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

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

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

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

2398:   Collective on sf

2400:   Input Parameters:
2401: + sf - The SF
2402: - rootSection - Section defined on root space

2404:   Output Parameters:
2405: + remoteOffsets - root offsets in leaf storage, or NULL
2406: - leafSection - Section defined on the leaf space

2408:   Level: advanced

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

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

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

2460:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

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

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

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

2522: /*@C
2523:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

2525:   Collective on sf

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

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

2535:   Level: developer

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

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

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

2570:   Collective on sf

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

2578:   Output Parameters:
2579: - sectionSF - The new SF

2581:   Note: Either rootSection or remoteOffsets can be specified

2583:   Level: advanced

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

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

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

2628:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2629:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2630:       PetscInt loff, dof, d;

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

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

2651:   Collective on section

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

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

2661:   Level: advanced

2663: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2664: @*/
2665: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2666: {

2670:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2671:   section->clObj     = obj;
2672:   PetscSectionDestroy(&section->clSection);
2673:   ISDestroy(&section->clPoints);
2674:   section->clSection = clSection;
2675:   section->clPoints  = clPoints;
2676:   return(0);
2677: }

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

2682:   Collective on section

2684:   Input Parameters:
2685: + section   - The PetscSection
2686: - obj       - A PetscObject which serves as the key for this index

2688:   Output Parameters:
2689: + clSection - Section giving the size of the closure of each point
2690: - clPoints  - IS giving the points in each closure

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

2694:   Level: advanced

2696: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2697: @*/
2698: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2699: {
2701:   if (section->clObj == obj) {
2702:     if (clSection) *clSection = section->clSection;
2703:     if (clPoints)  *clPoints  = section->clPoints;
2704:   } else {
2705:     if (clSection) *clSection = NULL;
2706:     if (clPoints)  *clPoints  = NULL;
2707:   }
2708:   return(0);
2709: }

2711: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2712: {
2713:   PetscInt       i;

2717:   if (section->clObj != obj) {
2718:     PetscSectionDestroy(&section->clSection);
2719:     ISDestroy(&section->clPoints);
2720:   }
2721:   section->clObj  = obj;
2722:   PetscFree(section->clPerm);
2723:   PetscFree(section->clInvPerm);
2724:   section->clSize = clSize;
2725:   if (mode == PETSC_COPY_VALUES) {
2726:     PetscMalloc1(clSize, &section->clPerm);
2727:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2728:     PetscArraycpy(section->clPerm, clPerm, clSize);
2729:   } else if (mode == PETSC_OWN_POINTER) {
2730:     section->clPerm = clPerm;
2731:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2732:   PetscMalloc1(clSize, &section->clInvPerm);
2733:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2734:   return(0);
2735: }

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

2740:   Not Collective

2742:   Input Parameters:
2743: + section - The PetscSection
2744: . obj     - A PetscObject which serves as the key for this index
2745: - perm    - Permutation of the cell dof closure

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

2750:   Level: intermediate

2752: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2753: @*/
2754: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2755: {
2756:   const PetscInt *clPerm = NULL;
2757:   PetscInt        clSize = 0;
2758:   PetscErrorCode  ierr;

2761:   if (perm) {
2762:     ISGetLocalSize(perm, &clSize);
2763:     ISGetIndices(perm, &clPerm);
2764:   }
2765:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2766:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2767:   return(0);
2768: }

2770: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2771: {
2773:   if (section->clObj == obj) {
2774:     if (size) *size = section->clSize;
2775:     if (perm) *perm = section->clPerm;
2776:   } else {
2777:     if (size) *size = 0;
2778:     if (perm) *perm = NULL;
2779:   }
2780:   return(0);
2781: }

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

2786:   Not collective

2788:   Input Parameters:
2789: + section   - The PetscSection
2790: - obj       - A PetscObject which serves as the key for this index

2792:   Output Parameter:
2793: . perm - The dof closure permutation

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

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

2800:   Level: intermediate

2802: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2803: @*/
2804: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2805: {
2806:   const PetscInt *clPerm;
2807:   PetscInt        clSize;
2808:   PetscErrorCode  ierr;

2811:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2812:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2813:   return(0);
2814: }

2816: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2817: {
2819:   if (section->clObj == obj) {
2820:     if (size) *size = section->clSize;
2821:     if (perm) *perm = section->clInvPerm;
2822:   } else {
2823:     if (size) *size = 0;
2824:     if (perm) *perm = NULL;
2825:   }
2826:   return(0);
2827: }

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

2832:   Not collective

2834:   Input Parameters:
2835: + section   - The PetscSection
2836: - obj       - A PetscObject which serves as the key for this index

2838:   Output Parameters:
2839: + size - The dof closure size
2840: - perm - The dof closure permutation

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

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

2847:   Level: intermediate

2849: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2850: @*/
2851: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2852: {
2853:   const PetscInt *clPerm;
2854:   PetscInt        clSize;
2855:   PetscErrorCode  ierr;

2858:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2859:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2860:   return(0);
2861: }

2863: /*@
2864:   PetscSectionGetField - Get the subsection associated with a single field

2866:   Input Parameters:
2867: + s     - The PetscSection
2868: - field - The field number

2870:   Output Parameter:
2871: . subs  - The subsection for the given field

2873:   Level: intermediate

2875: .seealso: PetscSectionSetNumFields()
2876: @*/
2877: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2878: {
2882:   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);
2883:   *subs = s->field[field];
2884:   return(0);
2885: }

2887: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2888: PetscFunctionList PetscSectionSymList = NULL;

2890: /*@
2891:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2893:   Collective

2895:   Input Parameter:
2896: . comm - the MPI communicator

2898:   Output Parameter:
2899: . sym - pointer to the new set of symmetries

2901:   Level: developer

2903: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2904: @*/
2905: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2906: {

2911:   ISInitializePackage();
2912:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2913:   return(0);
2914: }

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

2919:   Collective on PetscSectionSym

2921:   Input Parameters:
2922: + sym    - The section symmetry object
2923: - method - The name of the section symmetry type

2925:   Level: developer

2927: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2928: @*/
2929: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2930: {
2931:   PetscErrorCode (*r)(PetscSectionSym);
2932:   PetscBool      match;

2937:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2938:   if (match) return(0);

2940:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2941:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2942:   if (sym->ops->destroy) {
2943:     (*sym->ops->destroy)(sym);
2944:     sym->ops->destroy = NULL;
2945:   }
2946:   (*r)(sym);
2947:   PetscObjectChangeTypeName((PetscObject)sym,method);
2948:   return(0);
2949: }


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

2955:   Not Collective

2957:   Input Parameter:
2958: . sym  - The section symmetry

2960:   Output Parameter:
2961: . type - The index set type name

2963:   Level: developer

2965: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2966: @*/
2967: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2968: {
2972:   *type = ((PetscObject)sym)->type_name;
2973:   return(0);
2974: }

2976: /*@C
2977:   PetscSectionSymRegister - Adds a new section symmetry implementation

2979:   Not Collective

2981:   Input Parameters:
2982: + name        - The name of a new user-defined creation routine
2983: - create_func - The creation routine itself

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

2988:   Level: developer

2990: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2991: @*/
2992: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2993: {

2997:   ISInitializePackage();
2998:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2999:   return(0);
3000: }

3002: /*@
3003:    PetscSectionSymDestroy - Destroys a section symmetry.

3005:    Collective on PetscSectionSym

3007:    Input Parameters:
3008: .  sym - the section symmetry

3010:    Level: developer

3012: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3013: @*/
3014: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3015: {
3016:   SymWorkLink    link,next;

3020:   if (!*sym) return(0);
3022:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
3023:   if ((*sym)->ops->destroy) {
3024:     (*(*sym)->ops->destroy)(*sym);
3025:   }
3026:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3027:   for (link=(*sym)->workin; link; link=next) {
3028:     next = link->next;
3029:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3030:     PetscFree(link);
3031:   }
3032:   (*sym)->workin = NULL;
3033:   PetscHeaderDestroy(sym);
3034:   return(0);
3035: }

3037: /*@C
3038:    PetscSectionSymView - Displays a section symmetry

3040:    Collective on PetscSectionSym

3042:    Input Parameters:
3043: +  sym - the index set
3044: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

3046:    Level: developer

3048: .seealso: PetscViewerASCIIOpen()
3049: @*/
3050: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3051: {

3056:   if (!viewer) {
3057:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3058:   }
3061:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3062:   if (sym->ops->view) {
3063:     (*sym->ops->view)(sym,viewer);
3064:   }
3065:   return(0);
3066: }

3068: /*@
3069:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3071:   Collective on PetscSection

3073:   Input Parameters:
3074: + section - the section describing data layout
3075: - sym - the symmetry describing the affect of orientation on the access of the data

3077:   Level: developer

3079: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3080: @*/
3081: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3082: {

3087:   PetscSectionSymDestroy(&(section->sym));
3088:   if (sym) {
3091:     PetscObjectReference((PetscObject) sym);
3092:   }
3093:   section->sym = sym;
3094:   return(0);
3095: }

3097: /*@
3098:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3100:   Not collective

3102:   Input Parameters:
3103: . section - the section describing data layout

3105:   Output Parameters:
3106: . sym - the symmetry describing the affect of orientation on the access of the data

3108:   Level: developer

3110: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3111: @*/
3112: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3113: {
3116:   *sym = section->sym;
3117:   return(0);
3118: }

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

3123:   Collective on PetscSection

3125:   Input Parameters:
3126: + section - the section describing data layout
3127: . field - the field number
3128: - sym - the symmetry describing the affect of orientation on the access of the data

3130:   Level: developer

3132: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3133: @*/
3134: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3135: {

3140:   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);
3141:   PetscSectionSetSym(section->field[field],sym);
3142:   return(0);
3143: }

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

3148:   Collective on PetscSection

3150:   Input Parameters:
3151: + section - the section describing data layout
3152: - field - the field number

3154:   Output Parameters:
3155: . sym - the symmetry describing the affect of orientation on the access of the data

3157:   Level: developer

3159: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3160: @*/
3161: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3162: {
3165:   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);
3166:   *sym = section->field[field]->sym;
3167:   return(0);
3168: }

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

3173:   Not collective

3175:   Input Parameters:
3176: + section - the section
3177: . numPoints - the number of points
3178: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3179:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3180:     context, see DMPlexGetConeOrientation()).

3182:   Output Parameter:
3183: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3184: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3185:     identity).

3187:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3188: .vb
3189:      const PetscInt    **perms;
3190:      const PetscScalar **rots;
3191:      PetscInt            lOffset;

3193:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3194:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3195:        PetscInt           point = points[2*i], dof, sOffset;
3196:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3197:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3199:        PetscSectionGetDof(section,point,&dof);
3200:        PetscSectionGetOffset(section,point,&sOffset);

3202:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3203:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3204:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3205:        lOffset += dof;
3206:      }
3207:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3208: .ve

3210:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3211: .vb
3212:      const PetscInt    **perms;
3213:      const PetscScalar **rots;
3214:      PetscInt            lOffset;

3216:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3217:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3218:        PetscInt           point = points[2*i], dof, sOffset;
3219:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3220:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3222:        PetscSectionGetDof(section,point,&dof);
3223:        PetscSectionGetOffset(section,point,&sOff);

3225:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3226:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3227:        offset += dof;
3228:      }
3229:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3230: .ve

3232:   Level: developer

3234: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3235: @*/
3236: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3237: {
3238:   PetscSectionSym sym;
3239:   PetscErrorCode  ierr;

3244:   if (perms) *perms = NULL;
3245:   if (rots)  *rots  = NULL;
3246:   sym = section->sym;
3247:   if (sym && (perms || rots)) {
3248:     SymWorkLink link;

3250:     if (sym->workin) {
3251:       link        = sym->workin;
3252:       sym->workin = sym->workin->next;
3253:     } else {
3254:       PetscNewLog(sym,&link);
3255:     }
3256:     if (numPoints > link->numPoints) {
3257:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3258:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3259:       link->numPoints = numPoints;
3260:     }
3261:     link->next   = sym->workout;
3262:     sym->workout = link;
3263:     PetscArrayzero((PetscInt**)link->perms,numPoints);
3264:     PetscArrayzero((PetscInt**)link->rots,numPoints);
3265:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3266:     if (perms) *perms = link->perms;
3267:     if (rots)  *rots  = link->rots;
3268:   }
3269:   return(0);
3270: }

3272: /*@C
3273:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3275:   Not collective

3277:   Input Parameters:
3278: + section - the section
3279: . numPoints - the number of points
3280: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3281:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3282:     context, see DMPlexGetConeOrientation()).

3284:   Output Parameter:
3285: + perms - The permutations for the given orientations: set to NULL at conclusion
3286: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3288:   Level: developer

3290: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3291: @*/
3292: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3293: {
3294:   PetscSectionSym sym;

3298:   sym = section->sym;
3299:   if (sym && (perms || rots)) {
3300:     SymWorkLink *p,link;

3302:     for (p=&sym->workout; (link=*p); p=&link->next) {
3303:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3304:         *p          = link->next;
3305:         link->next  = sym->workin;
3306:         sym->workin = link;
3307:         if (perms) *perms = NULL;
3308:         if (rots)  *rots  = NULL;
3309:         return(0);
3310:       }
3311:     }
3312:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3313:   }
3314:   return(0);
3315: }

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

3320:   Not collective

3322:   Input Parameters:
3323: + section - the section
3324: . field - the field of the section
3325: . numPoints - the number of points
3326: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3327:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3328:     context, see DMPlexGetConeOrientation()).

3330:   Output Parameter:
3331: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3332: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3333:     identity).

3335:   Level: developer

3337: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3338: @*/
3339: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3340: {

3345:   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);
3346:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3347:   return(0);
3348: }

3350: /*@C
3351:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3353:   Not collective

3355:   Input Parameters:
3356: + section - the section
3357: . field - the field number
3358: . numPoints - the number of points
3359: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3360:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3361:     context, see DMPlexGetConeOrientation()).

3363:   Output Parameter:
3364: + perms - The permutations for the given orientations: set to NULL at conclusion
3365: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3367:   Level: developer

3369: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3370: @*/
3371: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3372: {

3377:   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);
3378:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3379:   return(0);
3380: }

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

3385:   Not collective

3387:   Input Parameter:
3388: . s - the global PetscSection

3390:   Output Parameters:
3391: . flg - the flag

3393:   Level: developer

3395: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3396: @*/
3397: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3398: {
3401:   *flg = s->useFieldOff;
3402:   return(0);
3403: }

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

3408:   Not collective

3410:   Input Parameters:
3411: + s   - the global PetscSection
3412: - flg - the flag

3414:   Level: developer

3416: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3417: @*/
3418: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3419: {
3422:   s->useFieldOff = flg;
3423:   return(0);
3424: }

3426: #define PetscSectionExpandPoints_Loop(TYPE) \
3427: { \
3428:   PetscInt i, n, o0, o1, size; \
3429:   TYPE *a0 = (TYPE*)origArray, *a1; \
3430:   PetscSectionGetStorageSize(s, &size); \
3431:   PetscMalloc1(size, &a1); \
3432:   for (i=0; i<npoints; i++) { \
3433:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3434:     PetscSectionGetOffset(s, i, &o1); \
3435:     PetscSectionGetDof(s, i, &n); \
3436:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3437:   } \
3438:   *newArray = (void*)a1; \
3439: }

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

3444:   Not collective

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

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

3456:   Level: developer

3458: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3459: @*/
3460: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3461: {
3462:   PetscSection        s;
3463:   const PetscInt      *points_;
3464:   PetscInt            i, n, npoints, pStart, pEnd;
3465:   PetscMPIInt         unitsize;
3466:   PetscErrorCode      ierr;

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