Actual source code: vsectionis.c

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

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

  9: PetscClassId PETSC_SECTION_CLASSID;

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

 14:   Collective on MPI_Comm

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

 20:   Level: developer

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

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

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

 43:   ISInitializePackage();

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

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

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

 71:   Collective on MPI_Comm

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

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

 79:   Level: developer

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

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

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

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

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

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

155:   Collective on MPI_Comm

157:   Input Parameter:
158: . section - the PetscSection

160:   Output Parameter:
161: . newSection - the copy

163:   Level: developer

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

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

179: /*@
180:   PetscSectionCompare - Compares two sections

182:   Collective on PetscSection

184:   Input Parameterss:
185: + s1 - the first PetscSection
186: - s2 - the second PetscSection

188:   Output Parameter:
189: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise 

191:   Level: developer

193:   Notes:
194:   Field names are disregarded.

196: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
197: @*/
198: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
199: {
200:   PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
201:   const PetscInt *idx1, *idx2;
202:   IS perm1, perm2;
203:   PetscBool flg;
204:   PetscMPIInt mflg;

211:   flg = PETSC_FALSE;

213:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
214:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
215:     *congruent = PETSC_FALSE;
216:     return(0);
217:   }

219:   PetscSectionGetChart(s1, &pStart, &pEnd);
220:   PetscSectionGetChart(s2, &n1, &n2);
221:   if (pStart != n1 || pEnd != n2) goto not_congruent;

223:   PetscSectionGetPermutation(s1, &perm1);
224:   PetscSectionGetPermutation(s2, &perm2);
225:   if (perm1 && perm2) {
226:     ISEqual(perm1, perm2, congruent);
227:     if (!(*congruent)) goto not_congruent;
228:   } else if (perm1 != perm2) goto not_congruent;

230:   for (p = pStart; p < pEnd; ++p) {
231:     PetscSectionGetOffset(s1, p, &n1);
232:     PetscSectionGetOffset(s2, p, &n2);
233:     if (n1 != n2) goto not_congruent;

235:     PetscSectionGetDof(s1, p, &n1);
236:     PetscSectionGetDof(s2, p, &n2);
237:     if (n1 != n2) goto not_congruent;

239:     PetscSectionGetConstraintDof(s1, p, &ncdof);
240:     PetscSectionGetConstraintDof(s2, p, &n2);
241:     if (ncdof != n2) goto not_congruent;

243:     PetscSectionGetConstraintIndices(s1, p, &idx1);
244:     PetscSectionGetConstraintIndices(s2, p, &idx2);
245:     PetscMemcmp(idx1, idx2, ncdof*sizeof(PetscInt), congruent);
246:     if (!(*congruent)) goto not_congruent;
247:   }

249:   PetscSectionGetNumFields(s1, &nfields);
250:   PetscSectionGetNumFields(s2, &n2);
251:   if (nfields != n2) goto not_congruent;

253:   for (f = 0; f < nfields; ++f) {
254:     PetscSectionGetFieldComponents(s1, f, &n1);
255:     PetscSectionGetFieldComponents(s2, f, &n2);
256:     if (n1 != n2) goto not_congruent;

258:     for (p = pStart; p < pEnd; ++p) {
259:       PetscSectionGetFieldOffset(s1, p, f, &n1);
260:       PetscSectionGetFieldOffset(s2, p, f, &n2);
261:       if (n1 != n2) goto not_congruent;

263:       PetscSectionGetFieldDof(s1, p, f, &n1);
264:       PetscSectionGetFieldDof(s2, p, f, &n2);
265:       if (n1 != n2) goto not_congruent;

267:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
268:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
269:       if (nfcdof != n2) goto not_congruent;

271:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
272:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
273:       PetscMemcmp(idx1, idx2, nfcdof*sizeof(PetscInt), congruent);
274:       if (!(*congruent)) goto not_congruent;
275:     }
276:   }

278:   flg = PETSC_TRUE;
279: not_congruent:
280:   MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
281:   return(0);
282: }

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

287:   Not collective

289:   Input Parameter:
290: . s - the PetscSection

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

295:   Level: intermediate

297: .seealso: PetscSectionSetNumFields()
298: @*/
299: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
300: {
304:   *numFields = s->numFields;
305:   return(0);
306: }

308: /*@
309:   PetscSectionSetNumFields - Sets the number of fields.

311:   Not collective

313:   Input Parameters:
314: + s - the PetscSection
315: - numFields - the number of fields

317:   Level: intermediate

319: .seealso: PetscSectionGetNumFields()
320: @*/
321: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
322: {
323:   PetscInt       f;

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

331:   s->numFields = numFields;
332:   PetscMalloc1(s->numFields, &s->numFieldComponents);
333:   PetscMalloc1(s->numFields, &s->fieldNames);
334:   PetscMalloc1(s->numFields, &s->field);
335:   for (f = 0; f < s->numFields; ++f) {
336:     char name[64];

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

340:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
341:     PetscSNPrintf(name, 64, "Field_%D", f);
342:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
343:   }
344:   return(0);
345: }

347: /*@C
348:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

350:   Not Collective

352:   Input Parameters:
353: + s     - the PetscSection
354: - field - the field number

356:   Output Parameter:
357: . fieldName - the field name

359:   Level: developer

361: .seealso: PetscSectionSetFieldName()
362: @*/
363: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
364: {
368:   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);
369:   *fieldName = s->fieldNames[field];
370:   return(0);
371: }

373: /*@C
374:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

376:   Not Collective

378:   Input Parameters:
379: + s     - the PetscSection
380: . field - the field number
381: - fieldName - the field name

383:   Level: developer

385: .seealso: PetscSectionGetFieldName()
386: @*/
387: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
388: {

394:   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);
395:   PetscFree(s->fieldNames[field]);
396:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
397:   return(0);
398: }

400: /*@
401:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

403:   Not collective

405:   Input Parameters:
406: + s - the PetscSection
407: - field - the field number

409:   Output Parameter:
410: . numComp - the number of field components

412:   Level: intermediate

414: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
415: @*/
416: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
417: {
421:   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);
422:   *numComp = s->numFieldComponents[field];
423:   return(0);
424: }

426: /*@
427:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

429:   Not collective

431:   Input Parameters:
432: + s - the PetscSection
433: . field - the field number
434: - numComp - the number of field components

436:   Level: intermediate

438: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
439: @*/
440: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
441: {
444:   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);
445:   s->numFieldComponents[field] = numComp;
446:   return(0);
447: }

449: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
450: {

454:   if (!s->bc) {
455:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
456:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
457:   }
458:   return(0);
459: }

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

464:   Not collective

466:   Input Parameter:
467: . s - the PetscSection

469:   Output Parameters:
470: + pStart - the first point
471: - pEnd - one past the last point

473:   Level: intermediate

475: .seealso: PetscSectionSetChart(), PetscSectionCreate()
476: @*/
477: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
478: {
481:   if (pStart) *pStart = s->pStart;
482:   if (pEnd)   *pEnd   = s->pEnd;
483:   return(0);
484: }

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

489:   Not collective

491:   Input Parameters:
492: + s - the PetscSection
493: . pStart - the first point
494: - pEnd - one past the last point

496:   Level: intermediate

498: .seealso: PetscSectionGetChart(), PetscSectionCreate()
499: @*/
500: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
501: {
502:   PetscInt       f;

507:   /* Cannot Reset() because it destroys field information */
508:   s->setup = PETSC_FALSE;
509:   PetscSectionDestroy(&s->bc);
510:   PetscFree(s->bcIndices);
511:   PetscFree2(s->atlasDof, s->atlasOff);

513:   s->pStart = pStart;
514:   s->pEnd   = pEnd;
515:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
516:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
517:   for (f = 0; f < s->numFields; ++f) {
518:     PetscSectionSetChart(s->field[f], pStart, pEnd);
519:   }
520:   return(0);
521: }

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

526:   Not collective

528:   Input Parameter:
529: . s - the PetscSection

531:   Output Parameters:
532: . perm - The permutation as an IS

534:   Level: intermediate

536: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
537: @*/
538: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
539: {
543:   return(0);
544: }

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

549:   Not collective

551:   Input Parameters:
552: + s - the PetscSection
553: - perm - the permutation of points

555:   Level: intermediate

557: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
558: @*/
559: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
560: {

566:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
567:   if (s->perm != perm) {
568:     ISDestroy(&s->perm);
569:     if (perm) {
570:       s->perm = perm;
571:       PetscObjectReference((PetscObject) s->perm);
572:     }
573:   }
574:   return(0);
575: }

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

580:   Not collective

582:   Input Parameters:
583: + s - the PetscSection
584: - point - the point

586:   Output Parameter:
587: . numDof - the number of dof

589:   Level: intermediate

591: .seealso: PetscSectionSetDof(), PetscSectionCreate()
592: @*/
593: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
594: {
598: #if defined(PETSC_USE_DEBUG)
599:   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);
600: #endif
601:   *numDof = s->atlasDof[point - s->pStart];
602:   return(0);
603: }

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

608:   Not collective

610:   Input Parameters:
611: + s - the PetscSection
612: . point - the point
613: - numDof - the number of dof

615:   Level: intermediate

617: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
618: @*/
619: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
620: {
623:   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);
624:   s->atlasDof[point - s->pStart] = numDof;
625:   return(0);
626: }

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

631:   Not collective

633:   Input Parameters:
634: + s - the PetscSection
635: . point - the point
636: - numDof - the number of additional dof

638:   Level: intermediate

640: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
641: @*/
642: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
643: {
646:   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);
647:   s->atlasDof[point - s->pStart] += numDof;
648:   return(0);
649: }

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

654:   Not collective

656:   Input Parameters:
657: + s - the PetscSection
658: . point - the point
659: - field - the field

661:   Output Parameter:
662: . numDof - the number of dof

664:   Level: intermediate

666: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
667: @*/
668: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
669: {

674:   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);
675:   PetscSectionGetDof(s->field[field], point, numDof);
676:   return(0);
677: }

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

682:   Not collective

684:   Input Parameters:
685: + s - the PetscSection
686: . point - the point
687: . field - the field
688: - numDof - the number of dof

690:   Level: intermediate

692: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
693: @*/
694: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
695: {

700:   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);
701:   PetscSectionSetDof(s->field[field], point, numDof);
702:   return(0);
703: }

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

708:   Not collective

710:   Input Parameters:
711: + s - the PetscSection
712: . point - the point
713: . field - the field
714: - numDof - the number of dof

716:   Level: intermediate

718: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
719: @*/
720: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
721: {

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

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

734:   Not collective

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

740:   Output Parameter:
741: . numDof - the number of dof which are fixed by constraints

743:   Level: intermediate

745: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
746: @*/
747: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
748: {

754:   if (s->bc) {
755:     PetscSectionGetDof(s->bc, point, numDof);
756:   } else *numDof = 0;
757:   return(0);
758: }

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

763:   Not collective

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

770:   Level: intermediate

772: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
773: @*/
774: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
775: {

780:   if (numDof) {
781:     PetscSectionCheckConstraints_Static(s);
782:     PetscSectionSetDof(s->bc, point, numDof);
783:   }
784:   return(0);
785: }

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

790:   Not collective

792:   Input Parameters:
793: + s - the PetscSection
794: . point - the point
795: - numDof - the number of additional dof which are fixed by constraints

797:   Level: intermediate

799: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
800: @*/
801: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
802: {

807:   if (numDof) {
808:     PetscSectionCheckConstraints_Static(s);
809:     PetscSectionAddDof(s->bc, point, numDof);
810:   }
811:   return(0);
812: }

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

817:   Not collective

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

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

827:   Level: intermediate

829: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
830: @*/
831: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
832: {

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

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

846:   Not collective

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

854:   Level: intermediate

856: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
857: @*/
858: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
859: {

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

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

872:   Not collective

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

880:   Level: intermediate

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

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

895: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
896: {

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

904:     PetscSectionSetUp(s->bc);
905:     PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
906:   }
907:   return(0);
908: }

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

913:   Not collective

915:   Input Parameter:
916: . s - the PetscSection

918:   Level: intermediate

920: .seealso: PetscSectionCreate()
921: @*/
922: PetscErrorCode PetscSectionSetUp(PetscSection s)
923: {
924:   const PetscInt *pind   = NULL;
925:   PetscInt        offset = 0, p, f;
926:   PetscErrorCode  ierr;

930:   if (s->setup) return(0);
931:   s->setup = PETSC_TRUE;
932:   if (s->perm) {ISGetIndices(s->perm, &pind);}
933:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
934:     const PetscInt q = pind ? pind[p] : p;

936:     s->atlasOff[q] = offset;
937:     offset        += s->atlasDof[q];
938:     s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
939:   }
940:   PetscSectionSetUpBC(s);
941:   /* Assume that all fields have the same chart */
942:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
943:     const PetscInt q   = pind ? pind[p] : p;
944:     PetscInt       off = s->atlasOff[q];

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

949:       sf->atlasOff[q] = off;
950:       off += sf->atlasDof[q];
951:     }
952:   }
953:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
954:   for (f = 0; f < s->numFields; ++f) {
955:     PetscSectionSetUpBC(s->field[f]);
956:   }
957:   return(0);
958: }

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

963:   Not collective

965:   Input Parameters:
966: . s - the PetscSection

968:   Output Parameter:
969: . maxDof - the maximum dof

971:   Level: intermediate

973: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
974: @*/
975: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
976: {
980:   *maxDof = s->maxDof;
981:   return(0);
982: }

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

987:   Not collective

989:   Input Parameters:
990: + s - the PetscSection
991: - size - the allocated size

993:   Output Parameter:
994: . size - the size of an array which can hold all the dofs

996:   Level: intermediate

998: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
999: @*/
1000: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1001: {
1002:   PetscInt p, n = 0;

1007:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1008:   *size = n;
1009:   return(0);
1010: }

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

1015:   Not collective

1017:   Input Parameters:
1018: + s - the PetscSection
1019: - point - the point

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

1024:   Level: intermediate

1026: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1027: @*/
1028: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1029: {
1030:   PetscInt p, n = 0;

1035:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1036:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1037:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1038:   }
1039:   *size = n;
1040:   return(0);
1041: }

1043: /*@
1044:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1045:   the local section and an SF describing the section point overlap.

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

1053:   Output Parameter:
1054:   . gsection - The PetscSection for the global field layout

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

1058:   Level: developer

1060: .seealso: PetscSectionCreate()
1061: @*/
1062: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1063: {
1064:   PetscSection    gs;
1065:   const PetscInt *pind = NULL;
1066:   PetscInt       *recv = NULL, *neg = NULL;
1067:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
1068:   PetscErrorCode  ierr;

1076:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1077:   PetscSectionGetChart(s, &pStart, &pEnd);
1078:   PetscSectionSetChart(gs, pStart, pEnd);
1079:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1080:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1081:   /* Must allocate for all points visible to SF, which may be more than this section */
1082:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1083:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
1084:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1085:     PetscMemzero(neg,nroots*sizeof(PetscInt));
1086:   }
1087:   /* Mark all local points with negative dof */
1088:   for (p = pStart; p < pEnd; ++p) {
1089:     PetscSectionGetDof(s, p, &dof);
1090:     PetscSectionSetDof(gs, p, dof);
1091:     PetscSectionGetConstraintDof(s, p, &cdof);
1092:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1093:     if (neg) neg[p] = -(dof+1);
1094:   }
1095:   PetscSectionSetUpBC(gs);
1096:   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]));}
1097:   if (nroots >= 0) {
1098:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1099:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1100:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1101:     for (p = pStart; p < pEnd; ++p) {
1102:       if (recv[p] < 0) {
1103:         gs->atlasDof[p-pStart] = recv[p];
1104:         PetscSectionGetDof(s, p, &dof);
1105:         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);
1106:       }
1107:     }
1108:   }
1109:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1110:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1111:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1112:     const PetscInt q = pind ? pind[p] : p;

1114:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1115:     gs->atlasOff[q] = off;
1116:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1117:   }
1118:   if (!localOffsets) {
1119:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1120:     globalOff -= off;
1121:   }
1122:   for (p = pStart, off = 0; p < pEnd; ++p) {
1123:     gs->atlasOff[p-pStart] += globalOff;
1124:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1125:   }
1126:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1127:   /* Put in negative offsets for ghost points */
1128:   if (nroots >= 0) {
1129:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1130:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1131:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1132:     for (p = pStart; p < pEnd; ++p) {
1133:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1134:     }
1135:   }
1136:   PetscFree2(neg,recv);
1137:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1138:   *gsection = gs;
1139:   return(0);
1140: }

1142: /*@
1143:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1144:   the local section and an SF describing the section point overlap.

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

1153:   Output Parameter:
1154:   . gsection - The PetscSection for the global field layout

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

1158:   Level: developer

1160: .seealso: PetscSectionCreate()
1161: @*/
1162: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1163: {
1164:   const PetscInt *pind = NULL;
1165:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1166:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1167:   PetscErrorCode  ierr;

1173:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1174:   PetscSectionGetChart(s, &pStart, &pEnd);
1175:   PetscSectionSetChart(*gsection, pStart, pEnd);
1176:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1177:   if (nroots >= 0) {
1178:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1179:     PetscCalloc1(nroots, &neg);
1180:     if (nroots > pEnd-pStart) {
1181:       PetscCalloc1(nroots, &tmpOff);
1182:     } else {
1183:       tmpOff = &(*gsection)->atlasDof[-pStart];
1184:     }
1185:   }
1186:   /* Mark ghost points with negative dof */
1187:   for (p = pStart; p < pEnd; ++p) {
1188:     for (e = 0; e < numExcludes; ++e) {
1189:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1190:         PetscSectionSetDof(*gsection, p, 0);
1191:         break;
1192:       }
1193:     }
1194:     if (e < numExcludes) continue;
1195:     PetscSectionGetDof(s, p, &dof);
1196:     PetscSectionSetDof(*gsection, p, dof);
1197:     PetscSectionGetConstraintDof(s, p, &cdof);
1198:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1199:     if (neg) neg[p] = -(dof+1);
1200:   }
1201:   PetscSectionSetUpBC(*gsection);
1202:   if (nroots >= 0) {
1203:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1204:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1205:     if (nroots > pEnd - pStart) {
1206:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1207:     }
1208:   }
1209:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1210:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1211:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1212:     const PetscInt q = pind ? pind[p] : p;

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

1239: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1240: {
1241:   PetscInt       pStart, pEnd, p, localSize = 0;

1245:   PetscSectionGetChart(s, &pStart, &pEnd);
1246:   for (p = pStart; p < pEnd; ++p) {
1247:     PetscInt dof;

1249:     PetscSectionGetDof(s, p, &dof);
1250:     if (dof > 0) ++localSize;
1251:   }
1252:   PetscLayoutCreate(comm, layout);
1253:   PetscLayoutSetLocalSize(*layout, localSize);
1254:   PetscLayoutSetBlockSize(*layout, 1);
1255:   PetscLayoutSetUp(*layout);
1256:   return(0);
1257: }

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

1262:   Input Parameters:
1263: + comm - The MPI_Comm
1264: - s    - The PetscSection

1266:   Output Parameter:
1267: . layout - The layout for the section

1269:   Level: developer

1271: .seealso: PetscSectionCreate()
1272: @*/
1273: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1274: {
1275:   PetscInt       pStart, pEnd, p, localSize = 0;

1281:   PetscSectionGetChart(s, &pStart, &pEnd);
1282:   for (p = pStart; p < pEnd; ++p) {
1283:     PetscInt dof,cdof;

1285:     PetscSectionGetDof(s, p, &dof);
1286:     PetscSectionGetConstraintDof(s, p, &cdof);
1287:     if (dof-cdof > 0) localSize += dof-cdof;
1288:   }
1289:   PetscLayoutCreate(comm, layout);
1290:   PetscLayoutSetLocalSize(*layout, localSize);
1291:   PetscLayoutSetBlockSize(*layout, 1);
1292:   PetscLayoutSetUp(*layout);
1293:   return(0);
1294: }

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

1299:   Not collective

1301:   Input Parameters:
1302: + s - the PetscSection
1303: - point - the point

1305:   Output Parameter:
1306: . offset - the offset

1308:   Level: intermediate

1310: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1311: @*/
1312: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1313: {
1317: #if defined(PETSC_USE_DEBUG)
1318:   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);
1319: #endif
1320:   *offset = s->atlasOff[point - s->pStart];
1321:   return(0);
1322: }

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

1327:   Not collective

1329:   Input Parameters:
1330: + s - the PetscSection
1331: . point - the point
1332: - offset - the offset

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

1336:   Level: intermediate

1338: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1339: @*/
1340: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1341: {
1344:   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);
1345:   s->atlasOff[point - s->pStart] = offset;
1346:   return(0);
1347: }

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

1352:   Not collective

1354:   Input Parameters:
1355: + s - the PetscSection
1356: . point - the point
1357: - field - the field

1359:   Output Parameter:
1360: . offset - the offset

1362:   Level: intermediate

1364: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1365: @*/
1366: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1367: {

1373:   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);
1374:   PetscSectionGetOffset(s->field[field], point, offset);
1375:   return(0);
1376: }

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

1381:   Not collective

1383:   Input Parameters:
1384: + s - the PetscSection
1385: . point - the point
1386: . field - the field
1387: - offset - the offset

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

1391:   Level: intermediate

1393: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1394: @*/
1395: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1396: {

1401:   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);
1402:   PetscSectionSetOffset(s->field[field], point, offset);
1403:   return(0);
1404: }

1406: /* This gives the offset on a point of the field, ignoring constraints */
1407: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1408: {
1409:   PetscInt       off, foff;

1415:   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);
1416:   PetscSectionGetOffset(s, point, &off);
1417:   PetscSectionGetOffset(s->field[field], point, &foff);
1418:   *offset = foff - off;
1419:   return(0);
1420: }

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

1425:   Not collective

1427:   Input Parameter:
1428: . s - the PetscSection

1430:   Output Parameters:
1431: + start - the minimum offset
1432: - end   - one more than the maximum offset

1434:   Level: intermediate

1436: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1437: @*/
1438: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1439: {
1440:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1445:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1446:   PetscSectionGetChart(s, &pStart, &pEnd);
1447:   for (p = 0; p < pEnd-pStart; ++p) {
1448:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1450:     if (off >= 0) {
1451:       os = PetscMin(os, off);
1452:       oe = PetscMax(oe, off+dof);
1453:     }
1454:   }
1455:   if (start) *start = os;
1456:   if (end)   *end   = oe;
1457:   return(0);
1458: }

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

1466:   if (!numFields) return(0);
1470:   PetscSectionGetNumFields(s, &nF);
1471:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1472:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1473:   PetscSectionSetNumFields(*subs, numFields);
1474:   for (f = 0; f < numFields; ++f) {
1475:     const char *name   = NULL;
1476:     PetscInt   numComp = 0;

1478:     PetscSectionGetFieldName(s, fields[f], &name);
1479:     PetscSectionSetFieldName(*subs, f, name);
1480:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1481:     PetscSectionSetFieldComponents(*subs, f, numComp);
1482:   }
1483:   PetscSectionGetChart(s, &pStart, &pEnd);
1484:   PetscSectionSetChart(*subs, pStart, pEnd);
1485:   for (p = pStart; p < pEnd; ++p) {
1486:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1488:     for (f = 0; f < numFields; ++f) {
1489:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1490:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1491:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1492:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1493:       dof  += fdof;
1494:       cdof += cfdof;
1495:     }
1496:     PetscSectionSetDof(*subs, p, dof);
1497:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1498:     maxCdof = PetscMax(cdof, maxCdof);
1499:   }
1500:   PetscSectionSetUp(*subs);
1501:   if (maxCdof) {
1502:     PetscInt *indices;

1504:     PetscMalloc1(maxCdof, &indices);
1505:     for (p = pStart; p < pEnd; ++p) {
1506:       PetscInt cdof;

1508:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1509:       if (cdof) {
1510:         const PetscInt *oldIndices = NULL;
1511:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1516:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1517:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1518:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1519:           /* This can be sped up if we assume sorted fields */
1520:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1521:             PetscInt oldfdof = 0;
1522:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1523:             oldFoff += oldfdof;
1524:           }
1525:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1526:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1527:           numConst += cfdof;
1528:           fOff     += fdof;
1529:         }
1530:         PetscSectionSetConstraintIndices(*subs, p, indices);
1531:       }
1532:     }
1533:     PetscFree(indices);
1534:   }
1535:   return(0);
1536: }

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

1544:   if (!len) return(0);
1545:   for (i = 0; i < len; ++i) {
1546:     PetscInt nf, pStarti, pEndi;

1548:     PetscSectionGetNumFields(s[i], &nf);
1549:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1550:     pStart = PetscMin(pStart, pStarti);
1551:     pEnd   = PetscMax(pEnd,   pEndi);
1552:     Nf += nf;
1553:   }
1554:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1555:   PetscSectionSetNumFields(*supers, Nf);
1556:   for (i = 0, f = 0; i < len; ++i) {
1557:     PetscInt nf, fi;

1559:     PetscSectionGetNumFields(s[i], &nf);
1560:     for (fi = 0; fi < nf; ++fi, ++f) {
1561:       const char *name   = NULL;
1562:       PetscInt   numComp = 0;

1564:       PetscSectionGetFieldName(s[i], fi, &name);
1565:       PetscSectionSetFieldName(*supers, f, name);
1566:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1567:       PetscSectionSetFieldComponents(*supers, f, numComp);
1568:     }
1569:   }
1570:   PetscSectionSetChart(*supers, pStart, pEnd);
1571:   for (p = pStart; p < pEnd; ++p) {
1572:     PetscInt dof = 0, cdof = 0;

1574:     for (i = 0, f = 0; i < len; ++i) {
1575:       PetscInt nf, fi, pStarti, pEndi;
1576:       PetscInt fdof = 0, cfdof = 0;

1578:       PetscSectionGetNumFields(s[i], &nf);
1579:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1580:       if ((p < pStarti) || (p >= pEndi)) continue;
1581:       for (fi = 0; fi < nf; ++fi, ++f) {
1582:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1583:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1584:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1585:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1586:         dof  += fdof;
1587:         cdof += cfdof;
1588:       }
1589:     }
1590:     PetscSectionSetDof(*supers, p, dof);
1591:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1592:     maxCdof = PetscMax(cdof, maxCdof);
1593:   }
1594:   PetscSectionSetUp(*supers);
1595:   if (maxCdof) {
1596:     PetscInt *indices;

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

1602:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1603:       if (cdof) {
1604:         PetscInt dof, numConst = 0, fOff = 0;

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

1610:           PetscSectionGetNumFields(s[i], &nf);
1611:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1612:           if ((p < pStarti) || (p >= pEndi)) continue;
1613:           for (fi = 0; fi < nf; ++fi, ++f) {
1614:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1615:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1616:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1617:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1618:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1619:             numConst += cfdof;
1620:           }
1621:           PetscSectionGetDof(s[i], p, &dof);
1622:           fOff += dof;
1623:         }
1624:         PetscSectionSetConstraintIndices(*supers, p, indices);
1625:       }
1626:     }
1627:     PetscFree(indices);
1628:   }
1629:   return(0);
1630: }

1632: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1633: {
1634:   const PetscInt *points = NULL, *indices = NULL;
1635:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1642:   PetscSectionGetNumFields(s, &numFields);
1643:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1644:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1645:   for (f = 0; f < numFields; ++f) {
1646:     const char *name   = NULL;
1647:     PetscInt   numComp = 0;

1649:     PetscSectionGetFieldName(s, f, &name);
1650:     PetscSectionSetFieldName(*subs, f, name);
1651:     PetscSectionGetFieldComponents(s, f, &numComp);
1652:     PetscSectionSetFieldComponents(*subs, f, numComp);
1653:   }
1654:   /* For right now, we do not try to squeeze the subchart */
1655:   if (subpointMap) {
1656:     ISGetSize(subpointMap, &numSubpoints);
1657:     ISGetIndices(subpointMap, &points);
1658:   }
1659:   PetscSectionGetChart(s, &pStart, &pEnd);
1660:   PetscSectionSetChart(*subs, 0, numSubpoints);
1661:   for (p = pStart; p < pEnd; ++p) {
1662:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1664:     PetscFindInt(p, numSubpoints, points, &subp);
1665:     if (subp < 0) continue;
1666:     for (f = 0; f < numFields; ++f) {
1667:       PetscSectionGetFieldDof(s, p, f, &fdof);
1668:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1669:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1670:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1671:     }
1672:     PetscSectionGetDof(s, p, &dof);
1673:     PetscSectionSetDof(*subs, subp, dof);
1674:     PetscSectionGetConstraintDof(s, p, &cdof);
1675:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1676:   }
1677:   PetscSectionSetUp(*subs);
1678:   /* Change offsets to original offsets */
1679:   for (p = pStart; p < pEnd; ++p) {
1680:     PetscInt off, foff = 0;

1682:     PetscFindInt(p, numSubpoints, points, &subp);
1683:     if (subp < 0) continue;
1684:     for (f = 0; f < numFields; ++f) {
1685:       PetscSectionGetFieldOffset(s, p, f, &foff);
1686:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1687:     }
1688:     PetscSectionGetOffset(s, p, &off);
1689:     PetscSectionSetOffset(*subs, subp, off);
1690:   }
1691:   /* Copy constraint indices */
1692:   for (subp = 0; subp < numSubpoints; ++subp) {
1693:     PetscInt cdof;

1695:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1696:     if (cdof) {
1697:       for (f = 0; f < numFields; ++f) {
1698:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1699:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1700:       }
1701:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1702:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1703:     }
1704:   }
1705:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1706:   return(0);
1707: }

1709: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1710: {
1711:   PetscInt       p;
1712:   PetscMPIInt    rank;

1716:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1717:   PetscViewerASCIIPushSynchronized(viewer);
1718:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1719:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1720:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1721:       PetscInt b;

1723:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1724:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1725:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1726:       }
1727:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1728:     } else {
1729:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1730:     }
1731:   }
1732:   PetscViewerFlush(viewer);
1733:   PetscViewerASCIIPopSynchronized(viewer);
1734:   if (s->sym) {
1735:     PetscViewerASCIIPushTab(viewer);
1736:     PetscSectionSymView(s->sym,viewer);
1737:     PetscViewerASCIIPopTab(viewer);
1738:   }
1739:   return(0);
1740: }

1742: /*@C
1743:   PetscSectionView - Views a PetscSection

1745:   Collective on PetscSection

1747:   Input Parameters:
1748: + s - the PetscSection object to view
1749: - v - the viewer

1751:   Level: developer

1753: .seealso PetscSectionCreate(), PetscSectionDestroy()
1754: @*/
1755: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1756: {
1757:   PetscBool      isascii;
1758:   PetscInt       f;

1763:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1765:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1766:   if (isascii) {
1767:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1768:     if (s->numFields) {
1769:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1770:       for (f = 0; f < s->numFields; ++f) {
1771:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1772:         PetscSectionView_ASCII(s->field[f], viewer);
1773:       }
1774:     } else {
1775:       PetscSectionView_ASCII(s, viewer);
1776:     }
1777:   }
1778:   return(0);
1779: }

1781: /*@
1782:   PetscSectionReset - Frees all section data.

1784:   Not collective

1786:   Input Parameters:
1787: . s - the PetscSection

1789:   Level: developer

1791: .seealso: PetscSection, PetscSectionCreate()
1792: @*/
1793: PetscErrorCode PetscSectionReset(PetscSection s)
1794: {
1795:   PetscInt       f;

1800:   PetscFree(s->numFieldComponents);
1801:   for (f = 0; f < s->numFields; ++f) {
1802:     PetscSectionDestroy(&s->field[f]);
1803:     PetscFree(s->fieldNames[f]);
1804:   }
1805:   PetscFree(s->fieldNames);
1806:   PetscFree(s->field);
1807:   PetscSectionDestroy(&s->bc);
1808:   PetscFree(s->bcIndices);
1809:   PetscFree2(s->atlasDof, s->atlasOff);
1810:   PetscSectionDestroy(&s->clSection);
1811:   ISDestroy(&s->clPoints);
1812:   ISDestroy(&s->perm);
1813:   PetscFree(s->clPerm);
1814:   PetscFree(s->clInvPerm);
1815:   PetscSectionSymDestroy(&s->sym);

1817:   s->pStart    = -1;
1818:   s->pEnd      = -1;
1819:   s->maxDof    = 0;
1820:   s->setup     = PETSC_FALSE;
1821:   s->numFields = 0;
1822:   s->clObj     = NULL;
1823:   return(0);
1824: }

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

1829:   Not collective

1831:   Input Parameters:
1832: . s - the PetscSection

1834:   Level: developer

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

1839: .seealso: PetscSection, PetscSectionCreate()
1840: @*/
1841: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1842: {

1846:   if (!*s) return(0);
1848:   if (--((PetscObject)(*s))->refct > 0) {
1849:     *s = NULL;
1850:     return(0);
1851:   }
1852:   PetscSectionReset(*s);
1853:   PetscHeaderDestroy(s);
1854:   return(0);
1855: }

1857: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1858: {
1859:   const PetscInt p = point - s->pStart;

1863:   *values = &baseArray[s->atlasOff[p]];
1864:   return(0);
1865: }

1867: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1868: {
1869:   PetscInt       *array;
1870:   const PetscInt p           = point - s->pStart;
1871:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1872:   PetscInt       cDim        = 0;

1877:   PetscSectionGetConstraintDof(s, p, &cDim);
1878:   array = &baseArray[s->atlasOff[p]];
1879:   if (!cDim) {
1880:     if (orientation >= 0) {
1881:       const PetscInt dim = s->atlasDof[p];
1882:       PetscInt       i;

1884:       if (mode == INSERT_VALUES) {
1885:         for (i = 0; i < dim; ++i) array[i] = values[i];
1886:       } else {
1887:         for (i = 0; i < dim; ++i) array[i] += values[i];
1888:       }
1889:     } else {
1890:       PetscInt offset = 0;
1891:       PetscInt j      = -1, field, i;

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

1896:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1897:         offset += dim;
1898:       }
1899:     }
1900:   } else {
1901:     if (orientation >= 0) {
1902:       const PetscInt dim  = s->atlasDof[p];
1903:       PetscInt       cInd = 0, i;
1904:       const PetscInt *cDof;

1906:       PetscSectionGetConstraintIndices(s, point, &cDof);
1907:       if (mode == INSERT_VALUES) {
1908:         for (i = 0; i < dim; ++i) {
1909:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1910:           array[i] = values[i];
1911:         }
1912:       } else {
1913:         for (i = 0; i < dim; ++i) {
1914:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1915:           array[i] += values[i];
1916:         }
1917:       }
1918:     } else {
1919:       const PetscInt *cDof;
1920:       PetscInt       offset  = 0;
1921:       PetscInt       cOffset = 0;
1922:       PetscInt       j       = 0, field;

1924:       PetscSectionGetConstraintIndices(s, point, &cDof);
1925:       for (field = 0; field < s->numFields; ++field) {
1926:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1927:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1928:         const PetscInt sDim = dim - tDim;
1929:         PetscInt       cInd = 0, i,k;

1931:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1932:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1933:           array[j] = values[k];
1934:         }
1935:         offset  += dim;
1936:         cOffset += dim - tDim;
1937:       }
1938:     }
1939:   }
1940:   return(0);
1941: }

1943: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1944: {
1948:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1949:   return(0);
1950: }

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

1955:   Input Parameters:
1956: + s     - The PetscSection
1957: - point - The point

1959:   Output Parameter:
1960: . indices - The constrained dofs

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

1964:   Level: advanced

1966: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1967: @*/
1968: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1969: {

1974:   if (s->bc) {
1975:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1976:   } else *indices = NULL;
1977:   return(0);
1978: }

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

1983:   Input Parameters:
1984: + s     - The PetscSection
1985: . point - The point
1986: - indices - The constrained dofs

1988:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1990:   Level: advanced

1992: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1993: @*/
1994: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1995: {

2000:   if (s->bc) {
2001:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2002:   }
2003:   return(0);
2004: }

2006: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2007: {

2012:   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);
2013:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2014:   return(0);
2015: }

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

2023:   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);
2024:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2025:   return(0);
2026: }

2028: /*@
2029:   PetscSectionPermute - Reorder the section according to the input point permutation

2031:   Collective on PetscSection

2033:   Input Parameter:
2034: + section - The PetscSection object
2035: - perm - The point permutation, old point p becomes new point perm[p]

2037:   Output Parameter:
2038: . sectionNew - The permuted PetscSection

2040:   Level: intermediate

2042: .keywords: mesh
2043: .seealso: MatPermute()
2044: @*/
2045: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2046: {
2047:   PetscSection    s = section, sNew;
2048:   const PetscInt *perm;
2049:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
2050:   PetscErrorCode  ierr;

2056:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2057:   PetscSectionGetNumFields(s, &numFields);
2058:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2059:   for (f = 0; f < numFields; ++f) {
2060:     const char *name;
2061:     PetscInt    numComp;

2063:     PetscSectionGetFieldName(s, f, &name);
2064:     PetscSectionSetFieldName(sNew, f, name);
2065:     PetscSectionGetFieldComponents(s, f, &numComp);
2066:     PetscSectionSetFieldComponents(sNew, f, numComp);
2067:   }
2068:   ISGetLocalSize(permutation, &numPoints);
2069:   ISGetIndices(permutation, &perm);
2070:   PetscSectionGetChart(s, &pStart, &pEnd);
2071:   PetscSectionSetChart(sNew, pStart, pEnd);
2072:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
2073:   for (p = pStart; p < pEnd; ++p) {
2074:     PetscInt dof, cdof;

2076:     PetscSectionGetDof(s, p, &dof);
2077:     PetscSectionSetDof(sNew, perm[p], dof);
2078:     PetscSectionGetConstraintDof(s, p, &cdof);
2079:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2080:     for (f = 0; f < numFields; ++f) {
2081:       PetscSectionGetFieldDof(s, p, f, &dof);
2082:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2083:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2084:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2085:     }
2086:   }
2087:   PetscSectionSetUp(sNew);
2088:   for (p = pStart; p < pEnd; ++p) {
2089:     const PetscInt *cind;
2090:     PetscInt        cdof;

2092:     PetscSectionGetConstraintDof(s, p, &cdof);
2093:     if (cdof) {
2094:       PetscSectionGetConstraintIndices(s, p, &cind);
2095:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2096:     }
2097:     for (f = 0; f < numFields; ++f) {
2098:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2099:       if (cdof) {
2100:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2101:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2102:       }
2103:     }
2104:   }
2105:   ISRestoreIndices(permutation, &perm);
2106:   *sectionNew = sNew;
2107:   return(0);
2108: }

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

2113:   Collective

2115:   Input Parameters:
2116: + sf - The SF
2117: - rootSection - Section defined on root space

2119:   Output Parameters:
2120: + remoteOffsets - root offsets in leaf storage, or NULL
2121: - leafSection - Section defined on the leaf space

2123:   Level: intermediate

2125: .seealso: PetscSFCreate()
2126: @*/
2127: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2128: {
2129:   PetscSF        embedSF;
2130:   const PetscInt *ilocal, *indices;
2131:   IS             selected;
2132:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

2136:   PetscSectionGetNumFields(rootSection, &numFields);
2137:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2138:   for (f = 0; f < numFields; ++f) {
2139:     PetscInt numComp = 0;
2140:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2141:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2142:   }
2143:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2144:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2145:   rpEnd = PetscMin(rpEnd,nroots);
2146:   rpEnd = PetscMax(rpStart,rpEnd);
2147:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2148:   ISGetIndices(selected, &indices);
2149:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2150:   ISRestoreIndices(selected, &indices);
2151:   ISDestroy(&selected);
2152:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2153:   if (nleaves && ilocal) {
2154:     for (i = 0; i < nleaves; ++i) {
2155:       lpStart = PetscMin(lpStart, ilocal[i]);
2156:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
2157:     }
2158:     ++lpEnd;
2159:   } else {
2160:     lpStart = 0;
2161:     lpEnd   = nleaves;
2162:   }
2163:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
2164:   /* Could fuse these at the cost of a copy and extra allocation */
2165:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2166:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2167:   for (f = 0; f < numFields; ++f) {
2168:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2169:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2170:   }
2171:   if (remoteOffsets) {
2172:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2173:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2174:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2175:   }
2176:   PetscSFDestroy(&embedSF);
2177:   PetscSectionSetUp(leafSection);
2178:   return(0);
2179: }

2181: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2182: {
2183:   PetscSF         embedSF;
2184:   const PetscInt *indices;
2185:   IS              selected;
2186:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2187:   PetscErrorCode  ierr;

2190:   *remoteOffsets = NULL;
2191:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2192:   if (numRoots < 0) return(0);
2193:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2194:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2195:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2196:   ISGetIndices(selected, &indices);
2197:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2198:   ISRestoreIndices(selected, &indices);
2199:   ISDestroy(&selected);
2200:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2201:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2202:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2203:   PetscSFDestroy(&embedSF);
2204:   return(0);
2205: }

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

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

2216:   Output Parameters:
2217: - sectionSF - The new SF

2219:   Note: Either rootSection or remoteOffsets can be specified

2221:   Level: intermediate

2223: .seealso: PetscSFCreate()
2224: @*/
2225: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2226: {
2227:   MPI_Comm          comm;
2228:   const PetscInt    *localPoints;
2229:   const PetscSFNode *remotePoints;
2230:   PetscInt          lpStart, lpEnd;
2231:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2232:   PetscInt          *localIndices;
2233:   PetscSFNode       *remoteIndices;
2234:   PetscInt          i, ind;
2235:   PetscErrorCode    ierr;

2243:   PetscObjectGetComm((PetscObject)sf,&comm);
2244:   PetscSFCreate(comm, sectionSF);
2245:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2246:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2247:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2248:   if (numRoots < 0) return(0);
2249:   for (i = 0; i < numPoints; ++i) {
2250:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2251:     PetscInt dof;

2253:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2254:       PetscSectionGetDof(leafSection, localPoint, &dof);
2255:       numIndices += dof;
2256:     }
2257:   }
2258:   PetscMalloc1(numIndices, &localIndices);
2259:   PetscMalloc1(numIndices, &remoteIndices);
2260:   /* Create new index graph */
2261:   for (i = 0, ind = 0; i < numPoints; ++i) {
2262:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2263:     PetscInt rank       = remotePoints[i].rank;

2265:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2266:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2267:       PetscInt loff, dof, d;

2269:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2270:       PetscSectionGetDof(leafSection, localPoint, &dof);
2271:       for (d = 0; d < dof; ++d, ++ind) {
2272:         localIndices[ind]        = loff+d;
2273:         remoteIndices[ind].rank  = rank;
2274:         remoteIndices[ind].index = remoteOffset+d;
2275:       }
2276:     }
2277:   }
2278:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2279:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2280:   return(0);
2281: }

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

2286:   Input Parameters:
2287: + section   - The PetscSection
2288: . obj       - A PetscObject which serves as the key for this index
2289: . clSection - Section giving the size of the closure of each point
2290: - clPoints  - IS giving the points in each closure

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

2294:   Level: intermediate

2296: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2297: @*/
2298: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2299: {

2303:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2304:   section->clObj     = obj;
2305:   PetscSectionDestroy(&section->clSection);
2306:   ISDestroy(&section->clPoints);
2307:   section->clSection = clSection;
2308:   section->clPoints  = clPoints;
2309:   return(0);
2310: }

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

2315:   Input Parameters:
2316: + section   - The PetscSection
2317: - obj       - A PetscObject which serves as the key for this index

2319:   Output Parameters:
2320: + clSection - Section giving the size of the closure of each point
2321: - clPoints  - IS giving the points in each closure

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

2325:   Level: intermediate

2327: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2328: @*/
2329: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2330: {
2332:   if (section->clObj == obj) {
2333:     if (clSection) *clSection = section->clSection;
2334:     if (clPoints)  *clPoints  = section->clPoints;
2335:   } else {
2336:     if (clSection) *clSection = NULL;
2337:     if (clPoints)  *clPoints  = NULL;
2338:   }
2339:   return(0);
2340: }

2342: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2343: {
2344:   PetscInt       i;

2348:   if (section->clObj != obj) {
2349:     PetscSectionDestroy(&section->clSection);
2350:     ISDestroy(&section->clPoints);
2351:   }
2352:   section->clObj  = obj;
2353:   PetscFree(section->clPerm);
2354:   PetscFree(section->clInvPerm);
2355:   section->clSize = clSize;
2356:   if (mode == PETSC_COPY_VALUES) {
2357:     PetscMalloc1(clSize, &section->clPerm);
2358:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2359:     PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2360:   } else if (mode == PETSC_OWN_POINTER) {
2361:     section->clPerm = clPerm;
2362:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2363:   PetscMalloc1(clSize, &section->clInvPerm);
2364:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2365:   return(0);
2366: }

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

2371:   Input Parameters:
2372: + section - The PetscSection
2373: . obj     - A PetscObject which serves as the key for this index
2374: - perm    - Permutation of the cell dof closure

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

2379:   Level: intermediate

2381: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2382: @*/
2383: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2384: {
2385:   const PetscInt *clPerm = NULL;
2386:   PetscInt        clSize = 0;
2387:   PetscErrorCode  ierr;

2390:   if (perm) {
2391:     ISGetLocalSize(perm, &clSize);
2392:     ISGetIndices(perm, &clPerm);
2393:   }
2394:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2395:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2396:   return(0);
2397: }

2399: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2400: {
2402:   if (section->clObj == obj) {
2403:     if (size) *size = section->clSize;
2404:     if (perm) *perm = section->clPerm;
2405:   } else {
2406:     if (size) *size = 0;
2407:     if (perm) *perm = NULL;
2408:   }
2409:   return(0);
2410: }

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

2415:   Input Parameters:
2416: + section   - The PetscSection
2417: - obj       - A PetscObject which serves as the key for this index

2419:   Output Parameter:
2420: . perm - The dof closure permutation

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

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

2427:   Level: intermediate

2429: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2430: @*/
2431: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2432: {
2433:   const PetscInt *clPerm;
2434:   PetscInt        clSize;
2435:   PetscErrorCode  ierr;

2438:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2439:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2440:   return(0);
2441: }

2443: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2444: {
2446:   if (section->clObj == obj) {
2447:     if (size) *size = section->clSize;
2448:     if (perm) *perm = section->clInvPerm;
2449:   } else {
2450:     if (size) *size = 0;
2451:     if (perm) *perm = NULL;
2452:   }
2453:   return(0);
2454: }

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

2459:   Input Parameters:
2460: + section   - The PetscSection
2461: - obj       - A PetscObject which serves as the key for this index

2463:   Output Parameters:
2464: + size - The dof closure size
2465: - perm - The dof closure permutation

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

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

2472:   Level: intermediate

2474: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2475: @*/
2476: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2477: {
2478:   const PetscInt *clPerm;
2479:   PetscInt        clSize;
2480:   PetscErrorCode  ierr;

2483:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2484:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2485:   return(0);
2486: }

2488: /*@
2489:   PetscSectionGetField - Get the subsection associated with a single field

2491:   Input Parameters:
2492: + s     - The PetscSection
2493: - field - The field number

2495:   Output Parameter:
2496: . subs  - The subsection for the given field

2498:   Level: intermediate

2500: .seealso: PetscSectionSetNumFields()
2501: @*/
2502: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2503: {

2509:   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);
2510:   PetscObjectReference((PetscObject) s->field[field]);
2511:   *subs = s->field[field];
2512:   return(0);
2513: }

2515: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2516: PetscFunctionList PetscSectionSymList = NULL;

2518: /*@
2519:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2521:   Collective on MPI_Comm

2523:   Input Parameter:
2524: . comm - the MPI communicator

2526:   Output Parameter:
2527: . sym - pointer to the new set of symmetries

2529:   Level: developer

2531: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2532: @*/
2533: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2534: {

2539:   ISInitializePackage();
2540:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2541:   return(0);
2542: }

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

2547:   Collective on PetscSectionSym

2549:   Input Parameters:
2550: + sym    - The section symmetry object
2551: - method - The name of the section symmetry type

2553:   Level: developer

2555: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2556: @*/
2557: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2558: {
2559:   PetscErrorCode (*r)(PetscSectionSym);
2560:   PetscBool      match;

2565:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2566:   if (match) return(0);

2568:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2569:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2570:   if (sym->ops->destroy) {
2571:     (*sym->ops->destroy)(sym);
2572:     sym->ops->destroy = NULL;
2573:   }
2574:   (*r)(sym);
2575:   PetscObjectChangeTypeName((PetscObject)sym,method);
2576:   return(0);
2577: }


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

2583:   Not Collective

2585:   Input Parameter:
2586: . sym  - The section symmetry

2588:   Output Parameter:
2589: . type - The index set type name

2591:   Level: developer

2593: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2594: @*/
2595: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2596: {
2600:   *type = ((PetscObject)sym)->type_name;
2601:   return(0);
2602: }

2604: /*@C
2605:   PetscSectionSymRegister - Adds a new section symmetry implementation

2607:   Not Collective

2609:   Input Parameters:
2610: + name        - The name of a new user-defined creation routine
2611: - create_func - The creation routine itself

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

2616:   Level: developer

2618: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2619: @*/
2620: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2621: {

2625:   ISInitializePackage();
2626:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2627:   return(0);
2628: }

2630: /*@
2631:    PetscSectionSymDestroy - Destroys a section symmetry.

2633:    Collective on PetscSectionSym

2635:    Input Parameters:
2636: .  sym - the section symmetry

2638:    Level: developer

2640: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2641: @*/
2642: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2643: {
2644:   SymWorkLink    link,next;

2648:   if (!*sym) return(0);
2650:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2651:   if ((*sym)->ops->destroy) {
2652:     (*(*sym)->ops->destroy)(*sym);
2653:   }
2654:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2655:   for (link=(*sym)->workin; link; link=next) {
2656:     next = link->next;
2657:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2658:     PetscFree(link);
2659:   }
2660:   (*sym)->workin = NULL;
2661:   PetscHeaderDestroy(sym);
2662:   return(0);
2663: }

2665: /*@C
2666:    PetscSectionSymView - Displays a section symmetry

2668:    Collective on PetscSectionSym

2670:    Input Parameters:
2671: +  sym - the index set
2672: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2674:    Level: developer

2676: .seealso: PetscViewerASCIIOpen()
2677: @*/
2678: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2679: {

2684:   if (!viewer) {
2685:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2686:   }
2689:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2690:   if (sym->ops->view) {
2691:     (*sym->ops->view)(sym,viewer);
2692:   }
2693:   return(0);
2694: }

2696: /*@
2697:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2699:   Collective on PetscSection

2701:   Input Parameters:
2702: + section - the section describing data layout
2703: - sym - the symmetry describing the affect of orientation on the access of the data

2705:   Level: developer

2707: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2708: @*/
2709: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2710: {

2715:   PetscSectionSymDestroy(&(section->sym));
2716:   if (sym) {
2719:     PetscObjectReference((PetscObject) sym);
2720:   }
2721:   section->sym = sym;
2722:   return(0);
2723: }

2725: /*@
2726:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2728:   Not collective

2730:   Input Parameters:
2731: . section - the section describing data layout

2733:   Output Parameters:
2734: . sym - the symmetry describing the affect of orientation on the access of the data

2736:   Level: developer

2738: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2739: @*/
2740: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2741: {
2744:   *sym = section->sym;
2745:   return(0);
2746: }

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

2751:   Collective on PetscSection

2753:   Input Parameters:
2754: + section - the section describing data layout
2755: . field - the field number
2756: - sym - the symmetry describing the affect of orientation on the access of the data

2758:   Level: developer

2760: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2761: @*/
2762: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2763: {

2768:   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);
2769:   PetscSectionSetSym(section->field[field],sym);
2770:   return(0);
2771: }

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

2776:   Collective on PetscSection

2778:   Input Parameters:
2779: + section - the section describing data layout
2780: - field - the field number

2782:   Output Parameters:
2783: . sym - the symmetry describing the affect of orientation on the access of the data

2785:   Level: developer

2787: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2788: @*/
2789: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2790: {
2793:   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);
2794:   *sym = section->field[field]->sym;
2795:   return(0);
2796: }

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

2801:   Not collective

2803:   Input Parameters:
2804: + section - the section
2805: . numPoints - the number of points
2806: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2807:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2808:     context, see DMPlexGetConeOrientation()).

2810:   Output Parameter:
2811: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2812: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2813:     identity).

2815:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2816: .vb
2817:      const PetscInt    **perms;
2818:      const PetscScalar **rots;
2819:      PetscInt            lOffset;

2821:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2822:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2823:        PetscInt           point = points[2*i], dof, sOffset;
2824:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2825:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2827:        PetscSectionGetDof(section,point,&dof);
2828:        PetscSectionGetOffset(section,point,&sOffset);

2830:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2831:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2832:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2833:        lOffset += dof;
2834:      }
2835:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2836: .ve

2838:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2839: .vb
2840:      const PetscInt    **perms;
2841:      const PetscScalar **rots;
2842:      PetscInt            lOffset;

2844:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2845:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2846:        PetscInt           point = points[2*i], dof, sOffset;
2847:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2848:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2850:        PetscSectionGetDof(section,point,&dof);
2851:        PetscSectionGetOffset(section,point,&sOff);

2853:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2854:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2855:        offset += dof;
2856:      }
2857:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2858: .ve

2860:   Level: developer

2862: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2863: @*/
2864: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2865: {
2866:   PetscSectionSym sym;
2867:   PetscErrorCode  ierr;

2872:   if (perms) *perms = NULL;
2873:   if (rots)  *rots  = NULL;
2874:   sym = section->sym;
2875:   if (sym && (perms || rots)) {
2876:     SymWorkLink link;

2878:     if (sym->workin) {
2879:       link        = sym->workin;
2880:       sym->workin = sym->workin->next;
2881:     } else {
2882:       PetscNewLog(sym,&link);
2883:     }
2884:     if (numPoints > link->numPoints) {
2885:       PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2886:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2887:       link->numPoints = numPoints;
2888:     }
2889:     link->next   = sym->workout;
2890:     sym->workout = link;
2891:     PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2892:     PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2893:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2894:     if (perms) *perms = link->perms;
2895:     if (rots)  *rots  = link->rots;
2896:   }
2897:   return(0);
2898: }

2900: /*@C
2901:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

2903:   Not collective

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

2912:   Output Parameter:
2913: + perms - The permutations for the given orientations: set to NULL at conclusion
2914: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2916:   Level: developer

2918: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2919: @*/
2920: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2921: {
2922:   PetscSectionSym sym;

2926:   sym = section->sym;
2927:   if (sym && (perms || rots)) {
2928:     SymWorkLink *p,link;

2930:     for (p=&sym->workout; (link=*p); p=&link->next) {
2931:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2932:         *p          = link->next;
2933:         link->next  = sym->workin;
2934:         sym->workin = link;
2935:         if (perms) *perms = NULL;
2936:         if (rots)  *rots  = NULL;
2937:         return(0);
2938:       }
2939:     }
2940:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2941:   }
2942:   return(0);
2943: }

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

2948:   Not collective

2950:   Input Parameters:
2951: + section - the section
2952: . field - the field of the section
2953: . numPoints - the number of points
2954: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2955:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2956:     context, see DMPlexGetConeOrientation()).

2958:   Output Parameter:
2959: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2960: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2961:     identity).

2963:   Level: developer

2965: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2966: @*/
2967: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2968: {

2973:   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);
2974:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2975:   return(0);
2976: }

2978: /*@C
2979:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

2981:   Not collective

2983:   Input Parameters:
2984: + section - the section
2985: . field - the field number
2986: . numPoints - the number of points
2987: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2988:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2989:     context, see DMPlexGetConeOrientation()).

2991:   Output Parameter:
2992: + perms - The permutations for the given orientations: set to NULL at conclusion
2993: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2995:   Level: developer

2997: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2998: @*/
2999: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3000: {

3005:   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);
3006:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3007:   return(0);
3008: }