Actual source code: vsectionis.c

petsc-master 2018-02-22
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

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

  9: PetscClassId PETSC_SECTION_CLASSID;

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

 14:   Collective on MPI_Comm

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

 20:   Level: developer

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

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

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

 42:   ISInitializePackage();

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

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

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

 70:   Collective on MPI_Comm

 72:   Input Parameter:
 73: . section - the PetscSection

 75:   Output Parameter:
 76: . newSection - the copy

 78:   Level: developer

 80: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
 81: @*/
 82: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 83: {
 84:   IS             perm;
 85:   PetscInt       numFields, f, pStart, pEnd, p;

 91:   PetscSectionGetNumFields(section, &numFields);
 92:   if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
 93:   for (f = 0; f < numFields; ++f) {
 94:     const char *name   = NULL;
 95:     PetscInt   numComp = 0;

 97:     PetscSectionGetFieldName(section, f, &name);
 98:     PetscSectionSetFieldName(newSection, f, name);
 99:     PetscSectionGetFieldComponents(section, f, &numComp);
100:     PetscSectionSetFieldComponents(newSection, f, numComp);
101:   }
102:   PetscSectionGetChart(section, &pStart, &pEnd);
103:   PetscSectionSetChart(newSection, pStart, pEnd);
104:   PetscSectionGetPermutation(section, &perm);
105:   PetscSectionSetPermutation(newSection, perm);
106:   for (p = pStart; p < pEnd; ++p) {
107:     PetscInt dof, cdof, fcdof = 0;

109:     PetscSectionGetDof(section, p, &dof);
110:     PetscSectionSetDof(newSection, p, dof);
111:     PetscSectionGetConstraintDof(section, p, &cdof);
112:     if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
113:     for (f = 0; f < numFields; ++f) {
114:       PetscSectionGetFieldDof(section, p, f, &dof);
115:       PetscSectionSetFieldDof(newSection, p, f, dof);
116:       if (cdof) {
117:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
118:         if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
119:       }
120:     }
121:   }
122:   PetscSectionSetUp(newSection);
123:   for (p = pStart; p < pEnd; ++p) {
124:     PetscInt       off, cdof, fcdof = 0;
125:     const PetscInt *cInd;

127:     /* Must set offsets in case they do not agree with the prefix sums */
128:     PetscSectionGetOffset(section, p, &off);
129:     PetscSectionSetOffset(newSection, p, off);
130:     PetscSectionGetConstraintDof(section, p, &cdof);
131:     if (cdof) {
132:       PetscSectionGetConstraintIndices(section, p, &cInd);
133:       PetscSectionSetConstraintIndices(newSection, p, cInd);
134:       for (f = 0; f < numFields; ++f) {
135:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
136:         if (fcdof) {
137:           PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
138:           PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
139:         }
140:       }
141:     }
142:   }
143:   return(0);
144: }

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

149:   Collective on MPI_Comm

151:   Input Parameter:
152: . section - the PetscSection

154:   Output Parameter:
155: . newSection - the copy

157:   Level: developer

159: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
160: @*/
161: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
162: {

168:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
169:   PetscSectionCopy(section, *newSection);
170:   return(0);
171: }

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

176:   Not collective

178:   Input Parameter:
179: . s - the PetscSection

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

184:   Level: intermediate

186: .seealso: PetscSectionSetNumFields()
187: @*/
188: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
189: {
193:   *numFields = s->numFields;
194:   return(0);
195: }

197: /*@
198:   PetscSectionSetNumFields - Sets the number of fields.

200:   Not collective

202:   Input Parameters:
203: + s - the PetscSection
204: - numFields - the number of fields

206:   Level: intermediate

208: .seealso: PetscSectionGetNumFields()
209: @*/
210: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
211: {
212:   PetscInt       f;

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

220:   s->numFields = numFields;
221:   PetscMalloc1(s->numFields, &s->numFieldComponents);
222:   PetscMalloc1(s->numFields, &s->fieldNames);
223:   PetscMalloc1(s->numFields, &s->field);
224:   for (f = 0; f < s->numFields; ++f) {
225:     char name[64];

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

229:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
230:     PetscSNPrintf(name, 64, "Field_%D", f);
231:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
232:   }
233:   return(0);
234: }

236: /*@C
237:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

239:   Not Collective

241:   Input Parameters:
242: + s     - the PetscSection
243: - field - the field number

245:   Output Parameter:
246: . fieldName - the field name

248:   Level: developer

250: .seealso: PetscSectionSetFieldName()
251: @*/
252: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
253: {
257:   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);
258:   *fieldName = s->fieldNames[field];
259:   return(0);
260: }

262: /*@C
263:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

265:   Not Collective

267:   Input Parameters:
268: + s     - the PetscSection
269: . field - the field number
270: - fieldName - the field name

272:   Level: developer

274: .seealso: PetscSectionGetFieldName()
275: @*/
276: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
277: {

283:   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);
284:   PetscFree(s->fieldNames[field]);
285:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
286:   return(0);
287: }

289: /*@
290:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

292:   Not collective

294:   Input Parameters:
295: + s - the PetscSection
296: - field - the field number

298:   Output Parameter:
299: . numComp - the number of field components

301:   Level: intermediate

303: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
304: @*/
305: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
306: {
310:   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);
311:   *numComp = s->numFieldComponents[field];
312:   return(0);
313: }

315: /*@
316:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

318:   Not collective

320:   Input Parameters:
321: + s - the PetscSection
322: . field - the field number
323: - numComp - the number of field components

325:   Level: intermediate

327: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
328: @*/
329: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
330: {
333:   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);
334:   s->numFieldComponents[field] = numComp;
335:   return(0);
336: }

338: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
339: {

343:   if (!s->bc) {
344:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
345:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
346:   }
347:   return(0);
348: }

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

353:   Not collective

355:   Input Parameter:
356: . s - the PetscSection

358:   Output Parameters:
359: + pStart - the first point
360: - pEnd - one past the last point

362:   Level: intermediate

364: .seealso: PetscSectionSetChart(), PetscSectionCreate()
365: @*/
366: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
367: {
370:   if (pStart) *pStart = s->pStart;
371:   if (pEnd)   *pEnd   = s->pEnd;
372:   return(0);
373: }

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

378:   Not collective

380:   Input Parameters:
381: + s - the PetscSection
382: . pStart - the first point
383: - pEnd - one past the last point

385:   Level: intermediate

387: .seealso: PetscSectionGetChart(), PetscSectionCreate()
388: @*/
389: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
390: {
391:   PetscInt       f;

396:   /* Cannot Reset() because it destroys field information */
397:   s->setup = PETSC_FALSE;
398:   PetscSectionDestroy(&s->bc);
399:   PetscFree(s->bcIndices);
400:   PetscFree2(s->atlasDof, s->atlasOff);

402:   s->pStart = pStart;
403:   s->pEnd   = pEnd;
404:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
405:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
406:   for (f = 0; f < s->numFields; ++f) {
407:     PetscSectionSetChart(s->field[f], pStart, pEnd);
408:   }
409:   return(0);
410: }

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

415:   Not collective

417:   Input Parameter:
418: . s - the PetscSection

420:   Output Parameters:
421: . perm - The permutation as an IS

423:   Level: intermediate

425: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
426: @*/
427: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
428: {
432:   return(0);
433: }

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

438:   Not collective

440:   Input Parameters:
441: + s - the PetscSection
442: - perm - the permutation of points

444:   Level: intermediate

446: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
447: @*/
448: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
449: {

455:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
456:   if (s->perm != perm) {
457:     ISDestroy(&s->perm);
458:     if (perm) {
459:       s->perm = perm;
460:       PetscObjectReference((PetscObject) s->perm);
461:     }
462:   }
463:   return(0);
464: }

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

469:   Not collective

471:   Input Parameters:
472: + s - the PetscSection
473: - point - the point

475:   Output Parameter:
476: . numDof - the number of dof

478:   Level: intermediate

480: .seealso: PetscSectionSetDof(), PetscSectionCreate()
481: @*/
482: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
483: {
487: #if defined(PETSC_USE_DEBUG)
488:   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);
489: #endif
490:   *numDof = s->atlasDof[point - s->pStart];
491:   return(0);
492: }

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

497:   Not collective

499:   Input Parameters:
500: + s - the PetscSection
501: . point - the point
502: - numDof - the number of dof

504:   Level: intermediate

506: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
507: @*/
508: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
509: {
512:   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);
513:   s->atlasDof[point - s->pStart] = numDof;
514:   return(0);
515: }

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

520:   Not collective

522:   Input Parameters:
523: + s - the PetscSection
524: . point - the point
525: - numDof - the number of additional dof

527:   Level: intermediate

529: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
530: @*/
531: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
532: {
535:   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);
536:   s->atlasDof[point - s->pStart] += numDof;
537:   return(0);
538: }

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

543:   Not collective

545:   Input Parameters:
546: + s - the PetscSection
547: . point - the point
548: - field - the field

550:   Output Parameter:
551: . numDof - the number of dof

553:   Level: intermediate

555: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
556: @*/
557: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
558: {

563:   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);
564:   PetscSectionGetDof(s->field[field], point, numDof);
565:   return(0);
566: }

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

571:   Not collective

573:   Input Parameters:
574: + s - the PetscSection
575: . point - the point
576: . field - the field
577: - numDof - the number of dof

579:   Level: intermediate

581: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
582: @*/
583: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
584: {

589:   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);
590:   PetscSectionSetDof(s->field[field], point, numDof);
591:   return(0);
592: }

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

597:   Not collective

599:   Input Parameters:
600: + s - the PetscSection
601: . point - the point
602: . field - the field
603: - numDof - the number of dof

605:   Level: intermediate

607: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
608: @*/
609: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
610: {

615:   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);
616:   PetscSectionAddDof(s->field[field], point, numDof);
617:   return(0);
618: }

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

623:   Not collective

625:   Input Parameters:
626: + s - the PetscSection
627: - point - the point

629:   Output Parameter:
630: . numDof - the number of dof which are fixed by constraints

632:   Level: intermediate

634: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
635: @*/
636: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
637: {

643:   if (s->bc) {
644:     PetscSectionGetDof(s->bc, point, numDof);
645:   } else *numDof = 0;
646:   return(0);
647: }

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

652:   Not collective

654:   Input Parameters:
655: + s - the PetscSection
656: . point - the point
657: - numDof - the number of dof which are fixed by constraints

659:   Level: intermediate

661: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
662: @*/
663: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
664: {

669:   if (numDof) {
670:     PetscSectionCheckConstraints_Static(s);
671:     PetscSectionSetDof(s->bc, point, numDof);
672:   }
673:   return(0);
674: }

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

679:   Not collective

681:   Input Parameters:
682: + s - the PetscSection
683: . point - the point
684: - numDof - the number of additional dof which are fixed by constraints

686:   Level: intermediate

688: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
689: @*/
690: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
691: {

696:   if (numDof) {
697:     PetscSectionCheckConstraints_Static(s);
698:     PetscSectionAddDof(s->bc, point, numDof);
699:   }
700:   return(0);
701: }

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

706:   Not collective

708:   Input Parameters:
709: + s - the PetscSection
710: . point - the point
711: - field - the field

713:   Output Parameter:
714: . numDof - the number of dof which are fixed by constraints

716:   Level: intermediate

718: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
719: @*/
720: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
721: {

727:   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);
728:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
729:   return(0);
730: }

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

735:   Not collective

737:   Input Parameters:
738: + s - the PetscSection
739: . point - the point
740: . field - the field
741: - numDof - the number of dof which are fixed by constraints

743:   Level: intermediate

745: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
746: @*/
747: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
748: {

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

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

761:   Not collective

763:   Input Parameters:
764: + s - the PetscSection
765: . point - the point
766: . field - the field
767: - numDof - the number of additional dof which are fixed by constraints

769:   Level: intermediate

771: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
772: @*/
773: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
774: {

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

784: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
785: {

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

793:     PetscSectionSetUp(s->bc);
794:     PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
795:   }
796:   return(0);
797: }

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

802:   Not collective

804:   Input Parameter:
805: . s - the PetscSection

807:   Level: intermediate

809: .seealso: PetscSectionCreate()
810: @*/
811: PetscErrorCode PetscSectionSetUp(PetscSection s)
812: {
813:   const PetscInt *pind   = NULL;
814:   PetscInt        offset = 0, p, f;
815:   PetscErrorCode  ierr;

819:   if (s->setup) return(0);
820:   s->setup = PETSC_TRUE;
821:   if (s->perm) {ISGetIndices(s->perm, &pind);}
822:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
823:     const PetscInt q = pind ? pind[p] : p;

825:     s->atlasOff[q] = offset;
826:     offset        += s->atlasDof[q];
827:     s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
828:   }
829:   PetscSectionSetUpBC(s);
830:   /* Assume that all fields have the same chart */
831:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
832:     const PetscInt q   = pind ? pind[p] : p;
833:     PetscInt       off = s->atlasOff[q];

835:     for (f = 0; f < s->numFields; ++f) {
836:       PetscSection sf = s->field[f];

838:       sf->atlasOff[q] = off;
839:       off += sf->atlasDof[q];
840:     }
841:   }
842:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
843:   for (f = 0; f < s->numFields; ++f) {
844:     PetscSectionSetUpBC(s->field[f]);
845:   }
846:   return(0);
847: }

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

852:   Not collective

854:   Input Parameters:
855: . s - the PetscSection

857:   Output Parameter:
858: . maxDof - the maximum dof

860:   Level: intermediate

862: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
863: @*/
864: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
865: {
869:   *maxDof = s->maxDof;
870:   return(0);
871: }

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

876:   Not collective

878:   Input Parameters:
879: + s - the PetscSection
880: - size - the allocated size

882:   Output Parameter:
883: . size - the size of an array which can hold all the dofs

885:   Level: intermediate

887: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
888: @*/
889: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
890: {
891:   PetscInt p, n = 0;

896:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
897:   *size = n;
898:   return(0);
899: }

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

904:   Not collective

906:   Input Parameters:
907: + s - the PetscSection
908: - point - the point

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

913:   Level: intermediate

915: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
916: @*/
917: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
918: {
919:   PetscInt p, n = 0;

924:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
925:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
926:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
927:   }
928:   *size = n;
929:   return(0);
930: }

932: /*@
933:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
934:   the local section and an SF describing the section point overlap.

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

942:   Output Parameter:
943:   . gsection - The PetscSection for the global field layout

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

947:   Level: developer

949: .seealso: PetscSectionCreate()
950: @*/
951: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
952: {
953:   PetscSection    gs;
954:   const PetscInt *pind = NULL;
955:   PetscInt       *recv = NULL, *neg = NULL;
956:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
957:   PetscErrorCode  ierr;

965:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
966:   PetscSectionGetChart(s, &pStart, &pEnd);
967:   PetscSectionSetChart(gs, pStart, pEnd);
968:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
969:   nlocal = nroots;              /* The local/leaf space matches global/root space */
970:   /* Must allocate for all points visible to SF, which may be more than this section */
971:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
972:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
973:     PetscMalloc2(nroots,&neg,nlocal,&recv);
974:     PetscMemzero(neg,nroots*sizeof(PetscInt));
975:   }
976:   /* Mark all local points with negative dof */
977:   for (p = pStart; p < pEnd; ++p) {
978:     PetscSectionGetDof(s, p, &dof);
979:     PetscSectionSetDof(gs, p, dof);
980:     PetscSectionGetConstraintDof(s, p, &cdof);
981:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
982:     if (neg) neg[p] = -(dof+1);
983:   }
984:   PetscSectionSetUpBC(gs);
985:   if (gs->bcIndices) {PetscMemcpy(gs->bcIndices, s->bcIndices, sizeof(PetscInt) * (gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]));}
986:   if (nroots >= 0) {
987:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
988:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
989:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
990:     for (p = pStart; p < pEnd; ++p) {
991:       if (recv[p] < 0) {
992:         gs->atlasDof[p-pStart] = recv[p];
993:         PetscSectionGetDof(s, p, &dof);
994:         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);
995:       }
996:     }
997:   }
998:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
999:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1000:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1001:     const PetscInt q = pind ? pind[p] : p;

1003:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1004:     gs->atlasOff[q] = off;
1005:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1006:   }
1007:   if (!localOffsets) {
1008:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1009:     globalOff -= off;
1010:   }
1011:   for (p = pStart, off = 0; p < pEnd; ++p) {
1012:     gs->atlasOff[p-pStart] += globalOff;
1013:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1014:   }
1015:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1016:   /* Put in negative offsets for ghost points */
1017:   if (nroots >= 0) {
1018:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1019:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1020:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1021:     for (p = pStart; p < pEnd; ++p) {
1022:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1023:     }
1024:   }
1025:   PetscFree2(neg,recv);
1026:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1027:   *gsection = gs;
1028:   return(0);
1029: }

1031: /*@
1032:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1033:   the local section and an SF describing the section point overlap.

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

1042:   Output Parameter:
1043:   . gsection - The PetscSection for the global field layout

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

1047:   Level: developer

1049: .seealso: PetscSectionCreate()
1050: @*/
1051: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1052: {
1053:   const PetscInt *pind = NULL;
1054:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1055:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1056:   PetscErrorCode  ierr;

1062:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1063:   PetscSectionGetChart(s, &pStart, &pEnd);
1064:   PetscSectionSetChart(*gsection, pStart, pEnd);
1065:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1066:   if (nroots >= 0) {
1067:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1068:     PetscCalloc1(nroots, &neg);
1069:     if (nroots > pEnd-pStart) {
1070:       PetscCalloc1(nroots, &tmpOff);
1071:     } else {
1072:       tmpOff = &(*gsection)->atlasDof[-pStart];
1073:     }
1074:   }
1075:   /* Mark ghost points with negative dof */
1076:   for (p = pStart; p < pEnd; ++p) {
1077:     for (e = 0; e < numExcludes; ++e) {
1078:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1079:         PetscSectionSetDof(*gsection, p, 0);
1080:         break;
1081:       }
1082:     }
1083:     if (e < numExcludes) continue;
1084:     PetscSectionGetDof(s, p, &dof);
1085:     PetscSectionSetDof(*gsection, p, dof);
1086:     PetscSectionGetConstraintDof(s, p, &cdof);
1087:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1088:     if (neg) neg[p] = -(dof+1);
1089:   }
1090:   PetscSectionSetUpBC(*gsection);
1091:   if (nroots >= 0) {
1092:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1093:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1094:     if (nroots > pEnd - pStart) {
1095:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1096:     }
1097:   }
1098:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1099:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1100:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1101:     const PetscInt q = pind ? pind[p] : p;

1103:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1104:     (*gsection)->atlasOff[q] = off;
1105:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1106:   }
1107:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1108:   globalOff -= off;
1109:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1110:     (*gsection)->atlasOff[p] += globalOff;
1111:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1112:   }
1113:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1114:   /* Put in negative offsets for ghost points */
1115:   if (nroots >= 0) {
1116:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1117:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1118:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1119:     if (nroots > pEnd - pStart) {
1120:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1121:     }
1122:   }
1123:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1124:   PetscFree(neg);
1125:   return(0);
1126: }

1128: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1129: {
1130:   PetscInt       pStart, pEnd, p, localSize = 0;

1134:   PetscSectionGetChart(s, &pStart, &pEnd);
1135:   for (p = pStart; p < pEnd; ++p) {
1136:     PetscInt dof;

1138:     PetscSectionGetDof(s, p, &dof);
1139:     if (dof > 0) ++localSize;
1140:   }
1141:   PetscLayoutCreate(comm, layout);
1142:   PetscLayoutSetLocalSize(*layout, localSize);
1143:   PetscLayoutSetBlockSize(*layout, 1);
1144:   PetscLayoutSetUp(*layout);
1145:   return(0);
1146: }

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

1151:   Input Parameters:
1152: + comm - The MPI_Comm
1153: - s    - The PetscSection

1155:   Output Parameter:
1156: . layout - The layout for the section

1158:   Level: developer

1160: .seealso: PetscSectionCreate()
1161: @*/
1162: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1163: {
1164:   PetscInt       pStart, pEnd, p, localSize = 0;

1170:   PetscSectionGetChart(s, &pStart, &pEnd);
1171:   for (p = pStart; p < pEnd; ++p) {
1172:     PetscInt dof,cdof;

1174:     PetscSectionGetDof(s, p, &dof);
1175:     PetscSectionGetConstraintDof(s, p, &cdof);
1176:     if (dof-cdof > 0) localSize += dof-cdof;
1177:   }
1178:   PetscLayoutCreate(comm, layout);
1179:   PetscLayoutSetLocalSize(*layout, localSize);
1180:   PetscLayoutSetBlockSize(*layout, 1);
1181:   PetscLayoutSetUp(*layout);
1182:   return(0);
1183: }

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

1188:   Not collective

1190:   Input Parameters:
1191: + s - the PetscSection
1192: - point - the point

1194:   Output Parameter:
1195: . offset - the offset

1197:   Level: intermediate

1199: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1200: @*/
1201: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1202: {
1206: #if defined(PETSC_USE_DEBUG)
1207:   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);
1208: #endif
1209:   *offset = s->atlasOff[point - s->pStart];
1210:   return(0);
1211: }

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

1216:   Not collective

1218:   Input Parameters:
1219: + s - the PetscSection
1220: . point - the point
1221: - offset - the offset

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

1225:   Level: intermediate

1227: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1228: @*/
1229: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1230: {
1233:   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);
1234:   s->atlasOff[point - s->pStart] = offset;
1235:   return(0);
1236: }

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

1241:   Not collective

1243:   Input Parameters:
1244: + s - the PetscSection
1245: . point - the point
1246: - field - the field

1248:   Output Parameter:
1249: . offset - the offset

1251:   Level: intermediate

1253: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1254: @*/
1255: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1256: {

1262:   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);
1263:   PetscSectionGetOffset(s->field[field], point, offset);
1264:   return(0);
1265: }

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

1270:   Not collective

1272:   Input Parameters:
1273: + s - the PetscSection
1274: . point - the point
1275: . field - the field
1276: - offset - the offset

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

1280:   Level: intermediate

1282: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1283: @*/
1284: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1285: {

1290:   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);
1291:   PetscSectionSetOffset(s->field[field], point, offset);
1292:   return(0);
1293: }

1295: /* This gives the offset on a point of the field, ignoring constraints */
1296: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1297: {
1298:   PetscInt       off, foff;

1304:   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);
1305:   PetscSectionGetOffset(s, point, &off);
1306:   PetscSectionGetOffset(s->field[field], point, &foff);
1307:   *offset = foff - off;
1308:   return(0);
1309: }

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

1314:   Not collective

1316:   Input Parameter:
1317: . s - the PetscSection

1319:   Output Parameters:
1320: + start - the minimum offset
1321: - end   - one more than the maximum offset

1323:   Level: intermediate

1325: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1326: @*/
1327: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1328: {
1329:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1334:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1335:   PetscSectionGetChart(s, &pStart, &pEnd);
1336:   for (p = 0; p < pEnd-pStart; ++p) {
1337:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1339:     if (off >= 0) {
1340:       os = PetscMin(os, off);
1341:       oe = PetscMax(oe, off+dof);
1342:     }
1343:   }
1344:   if (start) *start = os;
1345:   if (end)   *end   = oe;
1346:   return(0);
1347: }

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

1355:   if (!numFields) return(0);
1359:   PetscSectionGetNumFields(s, &nF);
1360:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1361:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1362:   PetscSectionSetNumFields(*subs, numFields);
1363:   for (f = 0; f < numFields; ++f) {
1364:     const char *name   = NULL;
1365:     PetscInt   numComp = 0;

1367:     PetscSectionGetFieldName(s, fields[f], &name);
1368:     PetscSectionSetFieldName(*subs, f, name);
1369:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1370:     PetscSectionSetFieldComponents(*subs, f, numComp);
1371:   }
1372:   PetscSectionGetChart(s, &pStart, &pEnd);
1373:   PetscSectionSetChart(*subs, pStart, pEnd);
1374:   for (p = pStart; p < pEnd; ++p) {
1375:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1377:     for (f = 0; f < numFields; ++f) {
1378:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1379:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1380:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1381:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1382:       dof  += fdof;
1383:       cdof += cfdof;
1384:     }
1385:     PetscSectionSetDof(*subs, p, dof);
1386:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1387:     maxCdof = PetscMax(cdof, maxCdof);
1388:   }
1389:   PetscSectionSetUp(*subs);
1390:   if (maxCdof) {
1391:     PetscInt *indices;

1393:     PetscMalloc1(maxCdof, &indices);
1394:     for (p = pStart; p < pEnd; ++p) {
1395:       PetscInt cdof;

1397:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1398:       if (cdof) {
1399:         const PetscInt *oldIndices = NULL;
1400:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1405:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1406:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1407:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1408:           /* This can be sped up if we assume sorted fields */
1409:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1410:             PetscInt oldfdof = 0;
1411:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1412:             oldFoff += oldfdof;
1413:           }
1414:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1415:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1416:           numConst += cfdof;
1417:           fOff     += fdof;
1418:         }
1419:         PetscSectionSetConstraintIndices(*subs, p, indices);
1420:       }
1421:     }
1422:     PetscFree(indices);
1423:   }
1424:   return(0);
1425: }

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

1433:   if (!len) return(0);
1434:   for (i = 0; i < len; ++i) {
1435:     PetscInt nf, pStarti, pEndi;

1437:     PetscSectionGetNumFields(s[i], &nf);
1438:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1439:     pStart = PetscMin(pStart, pStarti);
1440:     pEnd   = PetscMax(pEnd,   pEndi);
1441:     Nf += nf;
1442:   }
1443:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1444:   PetscSectionSetNumFields(*supers, Nf);
1445:   for (i = 0, f = 0; i < len; ++i) {
1446:     PetscInt nf, fi;

1448:     PetscSectionGetNumFields(s[i], &nf);
1449:     for (fi = 0; fi < nf; ++fi, ++f) {
1450:       const char *name   = NULL;
1451:       PetscInt   numComp = 0;

1453:       PetscSectionGetFieldName(s[i], fi, &name);
1454:       PetscSectionSetFieldName(*supers, f, name);
1455:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1456:       PetscSectionSetFieldComponents(*supers, f, numComp);
1457:     }
1458:   }
1459:   PetscSectionSetChart(*supers, pStart, pEnd);
1460:   for (p = pStart; p < pEnd; ++p) {
1461:     PetscInt dof = 0, cdof = 0;

1463:     for (i = 0, f = 0; i < len; ++i) {
1464:       PetscInt nf, fi, pStarti, pEndi;
1465:       PetscInt fdof = 0, cfdof = 0;

1467:       PetscSectionGetNumFields(s[i], &nf);
1468:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1469:       if ((p < pStarti) || (p >= pEndi)) continue;
1470:       for (fi = 0; fi < nf; ++fi, ++f) {
1471:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1472:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1473:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1474:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1475:         dof  += fdof;
1476:         cdof += cfdof;
1477:       }
1478:     }
1479:     PetscSectionSetDof(*supers, p, dof);
1480:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1481:     maxCdof = PetscMax(cdof, maxCdof);
1482:   }
1483:   PetscSectionSetUp(*supers);
1484:   if (maxCdof) {
1485:     PetscInt *indices;

1487:     PetscMalloc1(maxCdof, &indices);
1488:     for (p = pStart; p < pEnd; ++p) {
1489:       PetscInt cdof;

1491:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1492:       if (cdof) {
1493:         PetscInt dof, numConst = 0, fOff = 0;

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

1499:           PetscSectionGetNumFields(s[i], &nf);
1500:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1501:           if ((p < pStarti) || (p >= pEndi)) continue;
1502:           for (fi = 0; fi < nf; ++fi, ++f) {
1503:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1504:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1505:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1506:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1507:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1508:             numConst += cfdof;
1509:           }
1510:           PetscSectionGetDof(s[i], p, &dof);
1511:           fOff += dof;
1512:         }
1513:         PetscSectionSetConstraintIndices(*supers, p, indices);
1514:       }
1515:     }
1516:     PetscFree(indices);
1517:   }
1518:   return(0);
1519: }

1521: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1522: {
1523:   const PetscInt *points = NULL, *indices = NULL;
1524:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1531:   PetscSectionGetNumFields(s, &numFields);
1532:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1533:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1534:   for (f = 0; f < numFields; ++f) {
1535:     const char *name   = NULL;
1536:     PetscInt   numComp = 0;

1538:     PetscSectionGetFieldName(s, f, &name);
1539:     PetscSectionSetFieldName(*subs, f, name);
1540:     PetscSectionGetFieldComponents(s, f, &numComp);
1541:     PetscSectionSetFieldComponents(*subs, f, numComp);
1542:   }
1543:   /* For right now, we do not try to squeeze the subchart */
1544:   if (subpointMap) {
1545:     ISGetSize(subpointMap, &numSubpoints);
1546:     ISGetIndices(subpointMap, &points);
1547:   }
1548:   PetscSectionGetChart(s, &pStart, &pEnd);
1549:   PetscSectionSetChart(*subs, 0, numSubpoints);
1550:   for (p = pStart; p < pEnd; ++p) {
1551:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1553:     PetscFindInt(p, numSubpoints, points, &subp);
1554:     if (subp < 0) continue;
1555:     for (f = 0; f < numFields; ++f) {
1556:       PetscSectionGetFieldDof(s, p, f, &fdof);
1557:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1558:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1559:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1560:     }
1561:     PetscSectionGetDof(s, p, &dof);
1562:     PetscSectionSetDof(*subs, subp, dof);
1563:     PetscSectionGetConstraintDof(s, p, &cdof);
1564:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1565:   }
1566:   PetscSectionSetUp(*subs);
1567:   /* Change offsets to original offsets */
1568:   for (p = pStart; p < pEnd; ++p) {
1569:     PetscInt off, foff = 0;

1571:     PetscFindInt(p, numSubpoints, points, &subp);
1572:     if (subp < 0) continue;
1573:     for (f = 0; f < numFields; ++f) {
1574:       PetscSectionGetFieldOffset(s, p, f, &foff);
1575:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1576:     }
1577:     PetscSectionGetOffset(s, p, &off);
1578:     PetscSectionSetOffset(*subs, subp, off);
1579:   }
1580:   /* Copy constraint indices */
1581:   for (subp = 0; subp < numSubpoints; ++subp) {
1582:     PetscInt cdof;

1584:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1585:     if (cdof) {
1586:       for (f = 0; f < numFields; ++f) {
1587:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1588:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1589:       }
1590:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1591:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1592:     }
1593:   }
1594:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1595:   return(0);
1596: }

1598: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1599: {
1600:   PetscInt       p;
1601:   PetscMPIInt    rank;

1605:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1606:   PetscViewerASCIIPushSynchronized(viewer);
1607:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1608:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1609:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1610:       PetscInt b;

1612:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1613:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1614:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1615:       }
1616:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1617:     } else {
1618:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1619:     }
1620:   }
1621:   PetscViewerFlush(viewer);
1622:   PetscViewerASCIIPopSynchronized(viewer);
1623:   if (s->sym) {
1624:     PetscViewerASCIIPushTab(viewer);
1625:     PetscSectionSymView(s->sym,viewer);
1626:     PetscViewerASCIIPopTab(viewer);
1627:   }
1628:   return(0);
1629: }

1631: /*@C
1632:   PetscSectionView - Views a PetscSection

1634:   Collective on PetscSection

1636:   Input Parameters:
1637: + s - the PetscSection object to view
1638: - v - the viewer

1640:   Level: developer

1642: .seealso PetscSectionCreate(), PetscSectionDestroy()
1643: @*/
1644: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1645: {
1646:   PetscBool      isascii;
1647:   PetscInt       f;

1652:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1654:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1655:   if (isascii) {
1656:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1657:     if (s->numFields) {
1658:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1659:       for (f = 0; f < s->numFields; ++f) {
1660:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1661:         PetscSectionView_ASCII(s->field[f], viewer);
1662:       }
1663:     } else {
1664:       PetscSectionView_ASCII(s, viewer);
1665:     }
1666:   }
1667:   return(0);
1668: }

1670: /*@
1671:   PetscSectionReset - Frees all section data.

1673:   Not collective

1675:   Input Parameters:
1676: . s - the PetscSection

1678:   Level: developer

1680: .seealso: PetscSection, PetscSectionCreate()
1681: @*/
1682: PetscErrorCode PetscSectionReset(PetscSection s)
1683: {
1684:   PetscInt       f;

1689:   PetscFree(s->numFieldComponents);
1690:   for (f = 0; f < s->numFields; ++f) {
1691:     PetscSectionDestroy(&s->field[f]);
1692:     PetscFree(s->fieldNames[f]);
1693:   }
1694:   PetscFree(s->fieldNames);
1695:   PetscFree(s->field);
1696:   PetscSectionDestroy(&s->bc);
1697:   PetscFree(s->bcIndices);
1698:   PetscFree2(s->atlasDof, s->atlasOff);
1699:   PetscSectionDestroy(&s->clSection);
1700:   ISDestroy(&s->clPoints);
1701:   ISDestroy(&s->perm);
1702:   PetscFree(s->clPerm);
1703:   PetscFree(s->clInvPerm);
1704:   PetscSectionSymDestroy(&s->sym);

1706:   s->pStart    = -1;
1707:   s->pEnd      = -1;
1708:   s->maxDof    = 0;
1709:   s->setup     = PETSC_FALSE;
1710:   s->numFields = 0;
1711:   s->clObj     = NULL;
1712:   return(0);
1713: }

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

1718:   Not collective

1720:   Input Parameters:
1721: . s - the PetscSection

1723:   Level: developer

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

1728: .seealso: PetscSection, PetscSectionCreate()
1729: @*/
1730: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1731: {

1735:   if (!*s) return(0);
1737:   if (--((PetscObject)(*s))->refct > 0) {
1738:     *s = NULL;
1739:     return(0);
1740:   }
1741:   PetscSectionReset(*s);
1742:   PetscHeaderDestroy(s);
1743:   return(0);
1744: }

1746: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1747: {
1748:   const PetscInt p = point - s->pStart;

1752:   *values = &baseArray[s->atlasOff[p]];
1753:   return(0);
1754: }

1756: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1757: {
1758:   PetscInt       *array;
1759:   const PetscInt p           = point - s->pStart;
1760:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1761:   PetscInt       cDim        = 0;

1766:   PetscSectionGetConstraintDof(s, p, &cDim);
1767:   array = &baseArray[s->atlasOff[p]];
1768:   if (!cDim) {
1769:     if (orientation >= 0) {
1770:       const PetscInt dim = s->atlasDof[p];
1771:       PetscInt       i;

1773:       if (mode == INSERT_VALUES) {
1774:         for (i = 0; i < dim; ++i) array[i] = values[i];
1775:       } else {
1776:         for (i = 0; i < dim; ++i) array[i] += values[i];
1777:       }
1778:     } else {
1779:       PetscInt offset = 0;
1780:       PetscInt j      = -1, field, i;

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

1785:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1786:         offset += dim;
1787:       }
1788:     }
1789:   } else {
1790:     if (orientation >= 0) {
1791:       const PetscInt dim  = s->atlasDof[p];
1792:       PetscInt       cInd = 0, i;
1793:       const PetscInt *cDof;

1795:       PetscSectionGetConstraintIndices(s, point, &cDof);
1796:       if (mode == INSERT_VALUES) {
1797:         for (i = 0; i < dim; ++i) {
1798:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1799:           array[i] = values[i];
1800:         }
1801:       } else {
1802:         for (i = 0; i < dim; ++i) {
1803:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1804:           array[i] += values[i];
1805:         }
1806:       }
1807:     } else {
1808:       const PetscInt *cDof;
1809:       PetscInt       offset  = 0;
1810:       PetscInt       cOffset = 0;
1811:       PetscInt       j       = 0, field;

1813:       PetscSectionGetConstraintIndices(s, point, &cDof);
1814:       for (field = 0; field < s->numFields; ++field) {
1815:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1816:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1817:         const PetscInt sDim = dim - tDim;
1818:         PetscInt       cInd = 0, i,k;

1820:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1821:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1822:           array[j] = values[k];
1823:         }
1824:         offset  += dim;
1825:         cOffset += dim - tDim;
1826:       }
1827:     }
1828:   }
1829:   return(0);
1830: }

1832: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1833: {
1837:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1838:   return(0);
1839: }

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

1844:   Input Parameters:
1845: + s     - The PetscSection
1846: - point - The point

1848:   Output Parameter:
1849: . indices - The constrained dofs

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

1853:   Level: advanced

1855: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1856: @*/
1857: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1858: {

1863:   if (s->bc) {
1864:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1865:   } else *indices = NULL;
1866:   return(0);
1867: }

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

1872:   Input Parameters:
1873: + s     - The PetscSection
1874: . point - The point
1875: - indices - The constrained dofs

1877:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1879:   Level: advanced

1881: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1882: @*/
1883: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1884: {

1889:   if (s->bc) {
1890:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1891:   }
1892:   return(0);
1893: }

1895: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1896: {

1901:   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);
1902:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1903:   return(0);
1904: }

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

1912:   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);
1913:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1914:   return(0);
1915: }

1917: /*@
1918:   PetscSectionPermute - Reorder the section according to the input point permutation

1920:   Collective on PetscSection

1922:   Input Parameter:
1923: + section - The PetscSection object
1924: - perm - The point permutation, old point p becomes new point perm[p]

1926:   Output Parameter:
1927: . sectionNew - The permuted PetscSection

1929:   Level: intermediate

1931: .keywords: mesh
1932: .seealso: MatPermute()
1933: @*/
1934: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1935: {
1936:   PetscSection    s = section, sNew;
1937:   const PetscInt *perm;
1938:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1939:   PetscErrorCode  ierr;

1945:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1946:   PetscSectionGetNumFields(s, &numFields);
1947:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1948:   for (f = 0; f < numFields; ++f) {
1949:     const char *name;
1950:     PetscInt    numComp;

1952:     PetscSectionGetFieldName(s, f, &name);
1953:     PetscSectionSetFieldName(sNew, f, name);
1954:     PetscSectionGetFieldComponents(s, f, &numComp);
1955:     PetscSectionSetFieldComponents(sNew, f, numComp);
1956:   }
1957:   ISGetLocalSize(permutation, &numPoints);
1958:   ISGetIndices(permutation, &perm);
1959:   PetscSectionGetChart(s, &pStart, &pEnd);
1960:   PetscSectionSetChart(sNew, pStart, pEnd);
1961:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1962:   for (p = pStart; p < pEnd; ++p) {
1963:     PetscInt dof, cdof;

1965:     PetscSectionGetDof(s, p, &dof);
1966:     PetscSectionSetDof(sNew, perm[p], dof);
1967:     PetscSectionGetConstraintDof(s, p, &cdof);
1968:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1969:     for (f = 0; f < numFields; ++f) {
1970:       PetscSectionGetFieldDof(s, p, f, &dof);
1971:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1972:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1973:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1974:     }
1975:   }
1976:   PetscSectionSetUp(sNew);
1977:   for (p = pStart; p < pEnd; ++p) {
1978:     const PetscInt *cind;
1979:     PetscInt        cdof;

1981:     PetscSectionGetConstraintDof(s, p, &cdof);
1982:     if (cdof) {
1983:       PetscSectionGetConstraintIndices(s, p, &cind);
1984:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1985:     }
1986:     for (f = 0; f < numFields; ++f) {
1987:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1988:       if (cdof) {
1989:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1990:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1991:       }
1992:     }
1993:   }
1994:   ISRestoreIndices(permutation, &perm);
1995:   *sectionNew = sNew;
1996:   return(0);
1997: }

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

2002:   Collective

2004:   Input Parameters:
2005: + sf - The SF
2006: - rootSection - Section defined on root space

2008:   Output Parameters:
2009: + remoteOffsets - root offsets in leaf storage, or NULL
2010: - leafSection - Section defined on the leaf space

2012:   Level: intermediate

2014: .seealso: PetscSFCreate()
2015: @*/
2016: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2017: {
2018:   PetscSF        embedSF;
2019:   const PetscInt *ilocal, *indices;
2020:   IS             selected;
2021:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

2025:   PetscSectionGetNumFields(rootSection, &numFields);
2026:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2027:   for (f = 0; f < numFields; ++f) {
2028:     PetscInt numComp = 0;
2029:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2030:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2031:   }
2032:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2033:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2034:   rpEnd = PetscMin(rpEnd,nroots);
2035:   rpEnd = PetscMax(rpStart,rpEnd);
2036:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2037:   ISGetIndices(selected, &indices);
2038:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2039:   ISRestoreIndices(selected, &indices);
2040:   ISDestroy(&selected);
2041:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2042:   if (nleaves && ilocal) {
2043:     for (i = 0; i < nleaves; ++i) {
2044:       lpStart = PetscMin(lpStart, ilocal[i]);
2045:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
2046:     }
2047:     ++lpEnd;
2048:   } else {
2049:     lpStart = 0;
2050:     lpEnd   = nleaves;
2051:   }
2052:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
2053:   /* Could fuse these at the cost of a copy and extra allocation */
2054:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2055:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2056:   for (f = 0; f < numFields; ++f) {
2057:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2058:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2059:   }
2060:   if (remoteOffsets) {
2061:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2062:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2063:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2064:   }
2065:   PetscSFDestroy(&embedSF);
2066:   PetscSectionSetUp(leafSection);
2067:   return(0);
2068: }

2070: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2071: {
2072:   PetscSF         embedSF;
2073:   const PetscInt *indices;
2074:   IS              selected;
2075:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2076:   PetscErrorCode  ierr;

2079:   *remoteOffsets = NULL;
2080:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2081:   if (numRoots < 0) return(0);
2082:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2083:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2084:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2085:   ISGetIndices(selected, &indices);
2086:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2087:   ISRestoreIndices(selected, &indices);
2088:   ISDestroy(&selected);
2089:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2090:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2091:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2092:   PetscSFDestroy(&embedSF);
2093:   return(0);
2094: }

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

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

2105:   Output Parameters:
2106: - sectionSF - The new SF

2108:   Note: Either rootSection or remoteOffsets can be specified

2110:   Level: intermediate

2112: .seealso: PetscSFCreate()
2113: @*/
2114: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2115: {
2116:   MPI_Comm          comm;
2117:   const PetscInt    *localPoints;
2118:   const PetscSFNode *remotePoints;
2119:   PetscInt          lpStart, lpEnd;
2120:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2121:   PetscInt          *localIndices;
2122:   PetscSFNode       *remoteIndices;
2123:   PetscInt          i, ind;
2124:   PetscErrorCode    ierr;

2132:   PetscObjectGetComm((PetscObject)sf,&comm);
2133:   PetscSFCreate(comm, sectionSF);
2134:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2135:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2136:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2137:   if (numRoots < 0) return(0);
2138:   for (i = 0; i < numPoints; ++i) {
2139:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2140:     PetscInt dof;

2142:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2143:       PetscSectionGetDof(leafSection, localPoint, &dof);
2144:       numIndices += dof;
2145:     }
2146:   }
2147:   PetscMalloc1(numIndices, &localIndices);
2148:   PetscMalloc1(numIndices, &remoteIndices);
2149:   /* Create new index graph */
2150:   for (i = 0, ind = 0; i < numPoints; ++i) {
2151:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2152:     PetscInt rank       = remotePoints[i].rank;

2154:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2155:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2156:       PetscInt loff, dof, d;

2158:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2159:       PetscSectionGetDof(leafSection, localPoint, &dof);
2160:       for (d = 0; d < dof; ++d, ++ind) {
2161:         localIndices[ind]        = loff+d;
2162:         remoteIndices[ind].rank  = rank;
2163:         remoteIndices[ind].index = remoteOffset+d;
2164:       }
2165:     }
2166:   }
2167:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2168:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2169:   return(0);
2170: }

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

2175:   Input Parameters:
2176: + section   - The PetscSection
2177: . obj       - A PetscObject which serves as the key for this index
2178: . clSection - Section giving the size of the closure of each point
2179: - clPoints  - IS giving the points in each closure

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

2183:   Level: intermediate

2185: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2186: @*/
2187: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2188: {

2192:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2193:   section->clObj     = obj;
2194:   PetscSectionDestroy(&section->clSection);
2195:   ISDestroy(&section->clPoints);
2196:   section->clSection = clSection;
2197:   section->clPoints  = clPoints;
2198:   return(0);
2199: }

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

2204:   Input Parameters:
2205: + section   - The PetscSection
2206: - obj       - A PetscObject which serves as the key for this index

2208:   Output Parameters:
2209: + clSection - Section giving the size of the closure of each point
2210: - clPoints  - IS giving the points in each closure

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

2214:   Level: intermediate

2216: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2217: @*/
2218: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2219: {
2221:   if (section->clObj == obj) {
2222:     if (clSection) *clSection = section->clSection;
2223:     if (clPoints)  *clPoints  = section->clPoints;
2224:   } else {
2225:     if (clSection) *clSection = NULL;
2226:     if (clPoints)  *clPoints  = NULL;
2227:   }
2228:   return(0);
2229: }

2231: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2232: {
2233:   PetscInt       i;

2237:   if (section->clObj != obj) {
2238:     PetscSectionDestroy(&section->clSection);
2239:     ISDestroy(&section->clPoints);
2240:   }
2241:   section->clObj  = obj;
2242:   PetscFree(section->clPerm);
2243:   PetscFree(section->clInvPerm);
2244:   section->clSize = clSize;
2245:   if (mode == PETSC_COPY_VALUES) {
2246:     PetscMalloc1(clSize, &section->clPerm);
2247:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2248:     PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2249:   } else if (mode == PETSC_OWN_POINTER) {
2250:     section->clPerm = clPerm;
2251:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2252:   PetscMalloc1(clSize, &section->clInvPerm);
2253:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2254:   return(0);
2255: }

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

2260:   Input Parameters:
2261: + section - The PetscSection
2262: . obj     - A PetscObject which serves as the key for this index
2263: - perm    - Permutation of the cell dof closure

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

2268:   Level: intermediate

2270: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2271: @*/
2272: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2273: {
2274:   const PetscInt *clPerm = NULL;
2275:   PetscInt        clSize = 0;
2276:   PetscErrorCode  ierr;

2279:   if (perm) {
2280:     ISGetLocalSize(perm, &clSize);
2281:     ISGetIndices(perm, &clPerm);
2282:   }
2283:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2284:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2285:   return(0);
2286: }

2288: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2289: {
2291:   if (section->clObj == obj) {
2292:     if (size) *size = section->clSize;
2293:     if (perm) *perm = section->clPerm;
2294:   } else {
2295:     if (size) *size = 0;
2296:     if (perm) *perm = NULL;
2297:   }
2298:   return(0);
2299: }

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

2304:   Input Parameters:
2305: + section   - The PetscSection
2306: - obj       - A PetscObject which serves as the key for this index

2308:   Output Parameter:
2309: . perm - The dof closure permutation

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

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

2316:   Level: intermediate

2318: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2319: @*/
2320: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2321: {
2322:   const PetscInt *clPerm;
2323:   PetscInt        clSize;
2324:   PetscErrorCode  ierr;

2327:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2328:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2329:   return(0);
2330: }

2332: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2333: {
2335:   if (section->clObj == obj) {
2336:     if (size) *size = section->clSize;
2337:     if (perm) *perm = section->clInvPerm;
2338:   } else {
2339:     if (size) *size = 0;
2340:     if (perm) *perm = NULL;
2341:   }
2342:   return(0);
2343: }

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

2348:   Input Parameters:
2349: + section   - The PetscSection
2350: - obj       - A PetscObject which serves as the key for this index

2352:   Output Parameters:
2353: + size - The dof closure size
2354: - perm - The dof closure permutation

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

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

2361:   Level: intermediate

2363: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2364: @*/
2365: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2366: {
2367:   const PetscInt *clPerm;
2368:   PetscInt        clSize;
2369:   PetscErrorCode  ierr;

2372:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2373:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2374:   return(0);
2375: }

2377: /*@
2378:   PetscSectionGetField - Get the subsection associated with a single field

2380:   Input Parameters:
2381: + s     - The PetscSection
2382: - field - The field number

2384:   Output Parameter:
2385: . subs  - The subsection for the given field

2387:   Level: intermediate

2389: .seealso: PetscSectionSetNumFields()
2390: @*/
2391: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2392: {

2398:   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);
2399:   PetscObjectReference((PetscObject) s->field[field]);
2400:   *subs = s->field[field];
2401:   return(0);
2402: }

2404: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2405: PetscFunctionList PetscSectionSymList = NULL;

2407: /*@
2408:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2410:   Collective on MPI_Comm

2412:   Input Parameter:
2413: . comm - the MPI communicator

2415:   Output Parameter:
2416: . sym - pointer to the new set of symmetries

2418:   Level: developer

2420: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2421: @*/
2422: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2423: {

2428:   ISInitializePackage();
2429:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2430:   return(0);
2431: }

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

2436:   Collective on PetscSectionSym

2438:   Input Parameters:
2439: + sym    - The section symmetry object
2440: - method - The name of the section symmetry type

2442:   Level: developer

2444: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2445: @*/
2446: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2447: {
2448:   PetscErrorCode (*r)(PetscSectionSym);
2449:   PetscBool      match;

2454:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2455:   if (match) return(0);

2457:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2458:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2459:   if (sym->ops->destroy) {
2460:     (*sym->ops->destroy)(sym);
2461:     sym->ops->destroy = NULL;
2462:   }
2463:   (*r)(sym);
2464:   PetscObjectChangeTypeName((PetscObject)sym,method);
2465:   return(0);
2466: }


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

2472:   Not Collective

2474:   Input Parameter:
2475: . sym  - The section symmetry

2477:   Output Parameter:
2478: . type - The index set type name

2480:   Level: developer

2482: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2483: @*/
2484: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2485: {
2489:   *type = ((PetscObject)sym)->type_name;
2490:   return(0);
2491: }

2493: /*@C
2494:   PetscSectionSymRegister - Adds a new section symmetry implementation

2496:   Not Collective

2498:   Input Parameters:
2499: + name        - The name of a new user-defined creation routine
2500: - create_func - The creation routine itself

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

2505:   Level: developer

2507: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2508: @*/
2509: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2510: {

2514:   ISInitializePackage();
2515:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2516:   return(0);
2517: }

2519: /*@
2520:    PetscSectionSymDestroy - Destroys a section symmetry.

2522:    Collective on PetscSectionSym

2524:    Input Parameters:
2525: .  sym - the section symmetry

2527:    Level: developer

2529: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2530: @*/
2531: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2532: {
2533:   SymWorkLink    link,next;

2537:   if (!*sym) return(0);
2539:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2540:   if ((*sym)->ops->destroy) {
2541:     (*(*sym)->ops->destroy)(*sym);
2542:   }
2543:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2544:   for (link=(*sym)->workin; link; link=next) {
2545:     next = link->next;
2546:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2547:     PetscFree(link);
2548:   }
2549:   (*sym)->workin = NULL;
2550:   PetscHeaderDestroy(sym);
2551:   return(0);
2552: }

2554: /*@C
2555:    PetscSectionSymView - Displays a section symmetry

2557:    Collective on PetscSectionSym

2559:    Input Parameters:
2560: +  sym - the index set
2561: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2563:    Level: developer

2565: .seealso: PetscViewerASCIIOpen()
2566: @*/
2567: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2568: {

2573:   if (!viewer) {
2574:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2575:   }
2578:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2579:   if (sym->ops->view) {
2580:     (*sym->ops->view)(sym,viewer);
2581:   }
2582:   return(0);
2583: }

2585: /*@
2586:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2588:   Collective on PetscSection

2590:   Input Parameters:
2591: + section - the section describing data layout
2592: - sym - the symmetry describing the affect of orientation on the access of the data

2594:   Level: developer

2596: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2597: @*/
2598: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2599: {

2606:   PetscObjectReference((PetscObject)sym);
2607:   PetscSectionSymDestroy(&(section->sym));
2608:   section->sym = sym;
2609:   return(0);
2610: }

2612: /*@
2613:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2615:   Not collective

2617:   Input Parameters:
2618: . section - the section describing data layout

2620:   Output Parameters:
2621: . sym - the symmetry describing the affect of orientation on the access of the data

2623:   Level: developer

2625: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2626: @*/
2627: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2628: {
2631:   *sym = section->sym;
2632:   return(0);
2633: }

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

2638:   Collective on PetscSection

2640:   Input Parameters:
2641: + section - the section describing data layout
2642: . field - the field number
2643: - sym - the symmetry describing the affect of orientation on the access of the data

2645:   Level: developer

2647: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2648: @*/
2649: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2650: {

2655:   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);
2656:   PetscSectionSetSym(section->field[field],sym);
2657:   return(0);
2658: }

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

2663:   Collective on PetscSection

2665:   Input Parameters:
2666: + section - the section describing data layout
2667: - field - the field number

2669:   Output Parameters:
2670: . sym - the symmetry describing the affect of orientation on the access of the data

2672:   Level: developer

2674: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2675: @*/
2676: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2677: {
2680:   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);
2681:   *sym = section->field[field]->sym;
2682:   return(0);
2683: }

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

2688:   Not collective

2690:   Input Parameters:
2691: + section - the section
2692: . numPoints - the number of points
2693: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2694:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2695:     context, see DMPlexGetConeOrientation()).

2697:   Output Parameter:
2698: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2699: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2700:     identity).

2702:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2703: .vb
2704:      const PetscInt    **perms;
2705:      const PetscScalar **rots;
2706:      PetscInt            lOffset;

2708:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2709:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2710:        PetscInt           point = points[2*i], dof, sOffset;
2711:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2712:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2714:        PetscSectionGetDof(section,point,&dof);
2715:        PetscSectionGetOffset(section,point,&sOffset);

2717:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2718:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2719:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2720:        lOffset += dof;
2721:      }
2722:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2723: .ve

2725:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2726: .vb
2727:      const PetscInt    **perms;
2728:      const PetscScalar **rots;
2729:      PetscInt            lOffset;

2731:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2732:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2733:        PetscInt           point = points[2*i], dof, sOffset;
2734:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2735:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2737:        PetscSectionGetDof(section,point,&dof);
2738:        PetscSectionGetOffset(section,point,&sOff);

2740:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2741:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2742:        offset += dof;
2743:      }
2744:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2745: .ve

2747:   Level: developer

2749: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2750: @*/
2751: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2752: {
2753:   PetscSectionSym sym;
2754:   PetscErrorCode  ierr;

2759:   if (perms) *perms = NULL;
2760:   if (rots)  *rots  = NULL;
2761:   sym = section->sym;
2762:   if (sym && (perms || rots)) {
2763:     SymWorkLink link;

2765:     if (sym->workin) {
2766:       link        = sym->workin;
2767:       sym->workin = sym->workin->next;
2768:     } else {
2769:       PetscNewLog(sym,&link);
2770:     }
2771:     if (numPoints > link->numPoints) {
2772:       PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2773:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2774:       link->numPoints = numPoints;
2775:     }
2776:     link->next   = sym->workout;
2777:     sym->workout = link;
2778:     PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2779:     PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2780:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2781:     if (perms) *perms = link->perms;
2782:     if (rots)  *rots  = link->rots;
2783:   }
2784:   return(0);
2785: }

2787: /*@C
2788:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

2790:   Not collective

2792:   Input Parameters:
2793: + section - the section
2794: . numPoints - the number of points
2795: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2796:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2797:     context, see DMPlexGetConeOrientation()).

2799:   Output Parameter:
2800: + perms - The permutations for the given orientations: set to NULL at conclusion
2801: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2803:   Level: developer

2805: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2806: @*/
2807: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2808: {
2809:   PetscSectionSym sym;

2813:   sym = section->sym;
2814:   if (sym && (perms || rots)) {
2815:     SymWorkLink *p,link;

2817:     for (p=&sym->workout; (link=*p); p=&link->next) {
2818:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2819:         *p          = link->next;
2820:         link->next  = sym->workin;
2821:         sym->workin = link;
2822:         if (perms) *perms = NULL;
2823:         if (rots)  *rots  = NULL;
2824:         return(0);
2825:       }
2826:     }
2827:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2828:   }
2829:   return(0);
2830: }

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

2835:   Not collective

2837:   Input Parameters:
2838: + section - the section
2839: . field - the field of the section
2840: . numPoints - the number of points
2841: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2842:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2843:     context, see DMPlexGetConeOrientation()).

2845:   Output Parameter:
2846: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2847: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2848:     identity).

2850:   Level: developer

2852: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2853: @*/
2854: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2855: {

2860:   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);
2861:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2862:   return(0);
2863: }

2865: /*@C
2866:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

2868:   Not collective

2870:   Input Parameters:
2871: + section - the section
2872: . field - the field number
2873: . numPoints - the number of points
2874: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2875:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2876:     context, see DMPlexGetConeOrientation()).

2878:   Output Parameter:
2879: + perms - The permutations for the given orientations: set to NULL at conclusion
2880: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2882:   Level: developer

2884: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2885: @*/
2886: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2887: {

2892:   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);
2893:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
2894:   return(0);
2895: }