1 /*
2 This file contains routines for basic section object implementation.
3 */
4
5 #include <petsc/private/sectionimpl.h> /*I "petscsection.h" I*/
6 #include <petscsection.h>
7 #include <petscsf.h>
8 #include <petscviewer.h>
9
10 PetscClassId PETSC_SECTION_CLASSID;
11
12 /*@
13 PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
14
15 Collective
16
17 Input Parameters:
18 + comm - the MPI communicator
19 - s - pointer to the section
20
21 Level: beginner
22
23 Notes:
24 Typical calling sequence
25 $ PetscSectionCreate(MPI_Comm,PetscSection *);
26 $ PetscSectionSetNumFields(PetscSection, numFields);
27 $ PetscSectionSetChart(PetscSection,low,high);
28 $ PetscSectionSetDof(PetscSection,point,numdof);
29 $ PetscSectionSetUp(PetscSection);
30 $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
31 $ PetscSectionDestroy(PetscSection);
32
33 The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
34 recommended they not be used in user codes unless you really gain something in their use.
35
36 .seealso: PetscSection, PetscSectionDestroy()
37 @*/
PetscSectionCreate(MPI_Comm comm,PetscSection * s)38 PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
39 {
40 PetscErrorCode ierr;
41
42 PetscFunctionBegin;
43 PetscValidPointer(s,2);
44 ierr = ISInitializePackage();CHKERRQ(ierr);
45
46 ierr = PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);CHKERRQ(ierr);
47
48 (*s)->pStart = -1;
49 (*s)->pEnd = -1;
50 (*s)->perm = NULL;
51 (*s)->pointMajor = PETSC_TRUE;
52 (*s)->maxDof = 0;
53 (*s)->atlasDof = NULL;
54 (*s)->atlasOff = NULL;
55 (*s)->bc = NULL;
56 (*s)->bcIndices = NULL;
57 (*s)->setup = PETSC_FALSE;
58 (*s)->numFields = 0;
59 (*s)->fieldNames = NULL;
60 (*s)->field = NULL;
61 (*s)->useFieldOff = PETSC_FALSE;
62 (*s)->compNames = NULL;
63 (*s)->clObj = NULL;
64 (*s)->clHash = NULL;
65 (*s)->clSection = NULL;
66 (*s)->clPoints = NULL;
67 PetscFunctionReturn(0);
68 }
69
70 /*@
71 PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
72
73 Collective
74
75 Input Parameter:
76 . section - the PetscSection
77
78 Output Parameter:
79 . newSection - the copy
80
81 Level: intermediate
82
83 .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
84 @*/
PetscSectionCopy(PetscSection section,PetscSection newSection)85 PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86 {
87 PetscSectionSym sym;
88 IS perm;
89 PetscInt numFields, f, c, pStart, pEnd, p;
90 PetscErrorCode ierr;
91
92 PetscFunctionBegin;
93 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1);
94 PetscValidHeaderSpecific(newSection, PETSC_SECTION_CLASSID, 2);
95 ierr = PetscSectionReset(newSection);CHKERRQ(ierr);
96 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
97 if (numFields) {ierr = PetscSectionSetNumFields(newSection, numFields);CHKERRQ(ierr);}
98 for (f = 0; f < numFields; ++f) {
99 const char *fieldName = NULL, *compName = NULL;
100 PetscInt numComp = 0;
101
102 ierr = PetscSectionGetFieldName(section, f, &fieldName);CHKERRQ(ierr);
103 ierr = PetscSectionSetFieldName(newSection, f, fieldName);CHKERRQ(ierr);
104 ierr = PetscSectionGetFieldComponents(section, f, &numComp);CHKERRQ(ierr);
105 ierr = PetscSectionSetFieldComponents(newSection, f, numComp);CHKERRQ(ierr);
106 for (c = 0; c < numComp; ++c) {
107 ierr = PetscSectionGetComponentName(section, f, c, &compName);CHKERRQ(ierr);
108 ierr = PetscSectionSetComponentName(newSection, f, c, compName);CHKERRQ(ierr);
109 }
110 ierr = PetscSectionGetFieldSym(section, f, &sym);CHKERRQ(ierr);
111 ierr = PetscSectionSetFieldSym(newSection, f, sym);CHKERRQ(ierr);
112 }
113 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
114 ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
115 ierr = PetscSectionGetPermutation(section, &perm);CHKERRQ(ierr);
116 ierr = PetscSectionSetPermutation(newSection, perm);CHKERRQ(ierr);
117 ierr = PetscSectionGetSym(section, &sym);CHKERRQ(ierr);
118 ierr = PetscSectionSetSym(newSection, sym);CHKERRQ(ierr);
119 for (p = pStart; p < pEnd; ++p) {
120 PetscInt dof, cdof, fcdof = 0;
121
122 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
123 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
124 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
125 if (cdof) {ierr = PetscSectionSetConstraintDof(newSection, p, cdof);CHKERRQ(ierr);}
126 for (f = 0; f < numFields; ++f) {
127 ierr = PetscSectionGetFieldDof(section, p, f, &dof);CHKERRQ(ierr);
128 ierr = PetscSectionSetFieldDof(newSection, p, f, dof);CHKERRQ(ierr);
129 if (cdof) {
130 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
131 if (fcdof) {ierr = PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);CHKERRQ(ierr);}
132 }
133 }
134 }
135 ierr = PetscSectionSetUp(newSection);CHKERRQ(ierr);
136 for (p = pStart; p < pEnd; ++p) {
137 PetscInt off, cdof, fcdof = 0;
138 const PetscInt *cInd;
139
140 /* Must set offsets in case they do not agree with the prefix sums */
141 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
142 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
143 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
144 if (cdof) {
145 ierr = PetscSectionGetConstraintIndices(section, p, &cInd);CHKERRQ(ierr);
146 ierr = PetscSectionSetConstraintIndices(newSection, p, cInd);CHKERRQ(ierr);
147 for (f = 0; f < numFields; ++f) {
148 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
149 if (fcdof) {
150 ierr = PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);CHKERRQ(ierr);
151 ierr = PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);CHKERRQ(ierr);
152 }
153 }
154 }
155 }
156 PetscFunctionReturn(0);
157 }
158
159 /*@
160 PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
161
162 Collective
163
164 Input Parameter:
165 . section - the PetscSection
166
167 Output Parameter:
168 . newSection - the copy
169
170 Level: beginner
171
172 .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
173 @*/
PetscSectionClone(PetscSection section,PetscSection * newSection)174 PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
175 {
176 PetscErrorCode ierr;
177
178 PetscFunctionBegin;
179 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1);
180 PetscValidPointer(newSection, 2);
181 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);CHKERRQ(ierr);
182 ierr = PetscSectionCopy(section, *newSection);CHKERRQ(ierr);
183 PetscFunctionReturn(0);
184 }
185
186 /*@
187 PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
188
189 Collective on PetscSection
190
191 Input Parameter:
192 . section - the PetscSection
193
194 Options Database:
195 . -petscsection_point_major the dof order
196
197 Level: intermediate
198
199 .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
200 @*/
PetscSectionSetFromOptions(PetscSection s)201 PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
202 {
203 PetscErrorCode ierr;
204
205 PetscFunctionBegin;
206 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
207 ierr = PetscObjectOptionsBegin((PetscObject) s);CHKERRQ(ierr);
208 ierr = PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);CHKERRQ(ierr);
209 /* process any options handlers added with PetscObjectAddOptionsHandler() */
210 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);CHKERRQ(ierr);
211 ierr = PetscOptionsEnd();CHKERRQ(ierr);
212 ierr = PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");CHKERRQ(ierr);
213 PetscFunctionReturn(0);
214 }
215
216 /*@
217 PetscSectionCompare - Compares two sections
218
219 Collective on PetscSection
220
221 Input Parameters:
222 + s1 - the first PetscSection
223 - s2 - the second PetscSection
224
225 Output Parameter:
226 . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
227
228 Level: intermediate
229
230 Notes:
231 Field names are disregarded.
232
233 .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
234 @*/
PetscSectionCompare(PetscSection s1,PetscSection s2,PetscBool * congruent)235 PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
236 {
237 PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
238 const PetscInt *idx1, *idx2;
239 IS perm1, perm2;
240 PetscBool flg;
241 PetscMPIInt mflg;
242 PetscErrorCode ierr;
243
244 PetscFunctionBegin;
245 PetscValidHeaderSpecific(s1, PETSC_SECTION_CLASSID, 1);
246 PetscValidHeaderSpecific(s2, PETSC_SECTION_CLASSID, 2);
247 PetscValidIntPointer(congruent,3);
248 flg = PETSC_FALSE;
249
250 ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);CHKERRQ(ierr);
251 if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
252 *congruent = PETSC_FALSE;
253 PetscFunctionReturn(0);
254 }
255
256 ierr = PetscSectionGetChart(s1, &pStart, &pEnd);CHKERRQ(ierr);
257 ierr = PetscSectionGetChart(s2, &n1, &n2);CHKERRQ(ierr);
258 if (pStart != n1 || pEnd != n2) goto not_congruent;
259
260 ierr = PetscSectionGetPermutation(s1, &perm1);CHKERRQ(ierr);
261 ierr = PetscSectionGetPermutation(s2, &perm2);CHKERRQ(ierr);
262 if (perm1 && perm2) {
263 ierr = ISEqual(perm1, perm2, congruent);CHKERRQ(ierr);
264 if (!(*congruent)) goto not_congruent;
265 } else if (perm1 != perm2) goto not_congruent;
266
267 for (p = pStart; p < pEnd; ++p) {
268 ierr = PetscSectionGetOffset(s1, p, &n1);CHKERRQ(ierr);
269 ierr = PetscSectionGetOffset(s2, p, &n2);CHKERRQ(ierr);
270 if (n1 != n2) goto not_congruent;
271
272 ierr = PetscSectionGetDof(s1, p, &n1);CHKERRQ(ierr);
273 ierr = PetscSectionGetDof(s2, p, &n2);CHKERRQ(ierr);
274 if (n1 != n2) goto not_congruent;
275
276 ierr = PetscSectionGetConstraintDof(s1, p, &ncdof);CHKERRQ(ierr);
277 ierr = PetscSectionGetConstraintDof(s2, p, &n2);CHKERRQ(ierr);
278 if (ncdof != n2) goto not_congruent;
279
280 ierr = PetscSectionGetConstraintIndices(s1, p, &idx1);CHKERRQ(ierr);
281 ierr = PetscSectionGetConstraintIndices(s2, p, &idx2);CHKERRQ(ierr);
282 ierr = PetscArraycmp(idx1, idx2, ncdof, congruent);CHKERRQ(ierr);
283 if (!(*congruent)) goto not_congruent;
284 }
285
286 ierr = PetscSectionGetNumFields(s1, &nfields);CHKERRQ(ierr);
287 ierr = PetscSectionGetNumFields(s2, &n2);CHKERRQ(ierr);
288 if (nfields != n2) goto not_congruent;
289
290 for (f = 0; f < nfields; ++f) {
291 ierr = PetscSectionGetFieldComponents(s1, f, &n1);CHKERRQ(ierr);
292 ierr = PetscSectionGetFieldComponents(s2, f, &n2);CHKERRQ(ierr);
293 if (n1 != n2) goto not_congruent;
294
295 for (p = pStart; p < pEnd; ++p) {
296 ierr = PetscSectionGetFieldOffset(s1, p, f, &n1);CHKERRQ(ierr);
297 ierr = PetscSectionGetFieldOffset(s2, p, f, &n2);CHKERRQ(ierr);
298 if (n1 != n2) goto not_congruent;
299
300 ierr = PetscSectionGetFieldDof(s1, p, f, &n1);CHKERRQ(ierr);
301 ierr = PetscSectionGetFieldDof(s2, p, f, &n2);CHKERRQ(ierr);
302 if (n1 != n2) goto not_congruent;
303
304 ierr = PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);CHKERRQ(ierr);
305 ierr = PetscSectionGetFieldConstraintDof(s2, p, f, &n2);CHKERRQ(ierr);
306 if (nfcdof != n2) goto not_congruent;
307
308 ierr = PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);CHKERRQ(ierr);
309 ierr = PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);CHKERRQ(ierr);
310 ierr = PetscArraycmp(idx1, idx2, nfcdof, congruent);CHKERRQ(ierr);
311 if (!(*congruent)) goto not_congruent;
312 }
313 }
314
315 flg = PETSC_TRUE;
316 not_congruent:
317 ierr = MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));CHKERRQ(ierr);
318 PetscFunctionReturn(0);
319 }
320
321 /*@
322 PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
323
324 Not collective
325
326 Input Parameter:
327 . s - the PetscSection
328
329 Output Parameter:
330 . numFields - the number of fields defined, or 0 if none were defined
331
332 Level: intermediate
333
334 .seealso: PetscSectionSetNumFields()
335 @*/
PetscSectionGetNumFields(PetscSection s,PetscInt * numFields)336 PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
337 {
338 PetscFunctionBegin;
339 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
340 PetscValidPointer(numFields,2);
341 *numFields = s->numFields;
342 PetscFunctionReturn(0);
343 }
344
345 /*@
346 PetscSectionSetNumFields - Sets the number of fields.
347
348 Not collective
349
350 Input Parameters:
351 + s - the PetscSection
352 - numFields - the number of fields
353
354 Level: intermediate
355
356 .seealso: PetscSectionGetNumFields()
357 @*/
PetscSectionSetNumFields(PetscSection s,PetscInt numFields)358 PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
359 {
360 PetscInt f;
361 PetscErrorCode ierr;
362
363 PetscFunctionBegin;
364 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
365 if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
366 ierr = PetscSectionReset(s);CHKERRQ(ierr);
367
368 s->numFields = numFields;
369 ierr = PetscMalloc1(s->numFields, &s->numFieldComponents);CHKERRQ(ierr);
370 ierr = PetscMalloc1(s->numFields, &s->fieldNames);CHKERRQ(ierr);
371 ierr = PetscMalloc1(s->numFields, &s->compNames);CHKERRQ(ierr);
372 ierr = PetscMalloc1(s->numFields, &s->field);CHKERRQ(ierr);
373 for (f = 0; f < s->numFields; ++f) {
374 char name[64];
375
376 s->numFieldComponents[f] = 1;
377
378 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);CHKERRQ(ierr);
379 ierr = PetscSNPrintf(name, 64, "Field_%D", f);CHKERRQ(ierr);
380 ierr = PetscStrallocpy(name, (char **) &s->fieldNames[f]);CHKERRQ(ierr);
381 ierr = PetscSNPrintf(name, 64, "Component_0");CHKERRQ(ierr);
382 ierr = PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);CHKERRQ(ierr);
383 ierr = PetscStrallocpy(name, (char **) &s->compNames[f][0]);CHKERRQ(ierr);
384 }
385 PetscFunctionReturn(0);
386 }
387
388 /*@C
389 PetscSectionGetFieldName - Returns the name of a field in the PetscSection
390
391 Not Collective
392
393 Input Parameters:
394 + s - the PetscSection
395 - field - the field number
396
397 Output Parameter:
398 . fieldName - the field name
399
400 Level: intermediate
401
402 .seealso: PetscSectionSetFieldName()
403 @*/
PetscSectionGetFieldName(PetscSection s,PetscInt field,const char * fieldName[])404 PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
405 {
406 PetscFunctionBegin;
407 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
408 PetscValidPointer(fieldName, 3);
409 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);
410 *fieldName = s->fieldNames[field];
411 PetscFunctionReturn(0);
412 }
413
414 /*@C
415 PetscSectionSetFieldName - Sets the name of a field in the PetscSection
416
417 Not Collective
418
419 Input Parameters:
420 + s - the PetscSection
421 . field - the field number
422 - fieldName - the field name
423
424 Level: intermediate
425
426 .seealso: PetscSectionGetFieldName()
427 @*/
PetscSectionSetFieldName(PetscSection s,PetscInt field,const char fieldName[])428 PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
429 {
430 PetscErrorCode ierr;
431
432 PetscFunctionBegin;
433 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
434 if (fieldName) PetscValidCharPointer(fieldName, 3);
435 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);
436 ierr = PetscFree(s->fieldNames[field]);CHKERRQ(ierr);
437 ierr = PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);CHKERRQ(ierr);
438 PetscFunctionReturn(0);
439 }
440
441 /*@C
442 PetscSectionGetComponentName - Gets the name of a field component in the PetscSection
443
444 Not Collective
445
446 Input Parameters:
447 + s - the PetscSection
448 . field - the field number
449 . comp - the component number
450 - compName - the component name
451
452 Level: intermediate
453
454 .seealso: PetscSectionSetComponentName()
455 @*/
PetscSectionGetComponentName(PetscSection s,PetscInt field,PetscInt comp,const char * compName[])456 PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
457 {
458 PetscFunctionBegin;
459 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
460 PetscValidPointer(compName, 3);
461 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);
462 if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
463 *compName = s->compNames[field][comp];
464 PetscFunctionReturn(0);
465 }
466
467 /*@C
468 PetscSectionSetComponentName - Sets the name of a field component in the PetscSection
469
470 Not Collective
471
472 Input Parameters:
473 + s - the PetscSection
474 . field - the field number
475 . comp - the component number
476 - compName - the component name
477
478 Level: intermediate
479
480 .seealso: PetscSectionGetComponentName()
481 @*/
PetscSectionSetComponentName(PetscSection s,PetscInt field,PetscInt comp,const char compName[])482 PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
483 {
484 PetscErrorCode ierr;
485
486 PetscFunctionBegin;
487 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
488 if (compName) PetscValidCharPointer(compName, 3);
489 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);
490 if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
491 ierr = PetscFree(s->compNames[field][comp]);CHKERRQ(ierr);
492 ierr = PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);CHKERRQ(ierr);
493 PetscFunctionReturn(0);
494 }
495
496 /*@
497 PetscSectionGetFieldComponents - Returns the number of field components for the given field.
498
499 Not collective
500
501 Input Parameters:
502 + s - the PetscSection
503 - field - the field number
504
505 Output Parameter:
506 . numComp - the number of field components
507
508 Level: intermediate
509
510 .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
511 @*/
PetscSectionGetFieldComponents(PetscSection s,PetscInt field,PetscInt * numComp)512 PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
513 {
514 PetscFunctionBegin;
515 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
516 PetscValidPointer(numComp, 3);
517 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);
518 *numComp = s->numFieldComponents[field];
519 PetscFunctionReturn(0);
520 }
521
522 /*@
523 PetscSectionSetFieldComponents - Sets the number of field components for the given field.
524
525 Not collective
526
527 Input Parameters:
528 + s - the PetscSection
529 . field - the field number
530 - numComp - the number of field components
531
532 Level: intermediate
533
534 .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
535 @*/
PetscSectionSetFieldComponents(PetscSection s,PetscInt field,PetscInt numComp)536 PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
537 {
538 PetscErrorCode ierr;
539 PetscInt c;
540 char name[64];
541
542 PetscFunctionBegin;
543 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
544 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);
545 if (s->compNames) {
546 for (c = 0; c < s->numFieldComponents[field]; ++c) {
547 ierr = PetscFree(s->compNames[field][c]);CHKERRQ(ierr);
548 }
549 ierr = PetscFree(s->compNames[field]);CHKERRQ(ierr);
550 }
551
552 s->numFieldComponents[field] = numComp;
553 if (numComp) {
554 ierr = PetscMalloc1(numComp, (char ***) &s->compNames[field]);CHKERRQ(ierr);
555 for (c = 0; c < numComp; ++c) {
556 ierr = PetscSNPrintf(name, 64, "%D", c);CHKERRQ(ierr);
557 ierr = PetscStrallocpy(name, (char **) &s->compNames[field][c]);CHKERRQ(ierr);
558 }
559 }
560 PetscFunctionReturn(0);
561 }
562
PetscSectionCheckConstraints_Static(PetscSection s)563 static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
564 {
565 PetscErrorCode ierr;
566
567 PetscFunctionBegin;
568 if (!s->bc) {
569 ierr = PetscSectionCreate(PETSC_COMM_SELF, &s->bc);CHKERRQ(ierr);
570 ierr = PetscSectionSetChart(s->bc, s->pStart, s->pEnd);CHKERRQ(ierr);
571 }
572 PetscFunctionReturn(0);
573 }
574
575 /*@
576 PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
577
578 Not collective
579
580 Input Parameter:
581 . s - the PetscSection
582
583 Output Parameters:
584 + pStart - the first point
585 - pEnd - one past the last point
586
587 Level: intermediate
588
589 .seealso: PetscSectionSetChart(), PetscSectionCreate()
590 @*/
PetscSectionGetChart(PetscSection s,PetscInt * pStart,PetscInt * pEnd)591 PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
592 {
593 PetscFunctionBegin;
594 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
595 if (pStart) *pStart = s->pStart;
596 if (pEnd) *pEnd = s->pEnd;
597 PetscFunctionReturn(0);
598 }
599
600 /*@
601 PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
602
603 Not collective
604
605 Input Parameters:
606 + s - the PetscSection
607 . pStart - the first point
608 - pEnd - one past the last point
609
610 Level: intermediate
611
612 .seealso: PetscSectionGetChart(), PetscSectionCreate()
613 @*/
PetscSectionSetChart(PetscSection s,PetscInt pStart,PetscInt pEnd)614 PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
615 {
616 PetscInt f;
617 PetscErrorCode ierr;
618
619 PetscFunctionBegin;
620 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
621 /* Cannot Reset() because it destroys field information */
622 s->setup = PETSC_FALSE;
623 ierr = PetscSectionDestroy(&s->bc);CHKERRQ(ierr);
624 ierr = PetscFree(s->bcIndices);CHKERRQ(ierr);
625 ierr = PetscFree2(s->atlasDof, s->atlasOff);CHKERRQ(ierr);
626
627 s->pStart = pStart;
628 s->pEnd = pEnd;
629 ierr = PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);CHKERRQ(ierr);
630 ierr = PetscArrayzero(s->atlasDof, pEnd - pStart);CHKERRQ(ierr);
631 for (f = 0; f < s->numFields; ++f) {
632 ierr = PetscSectionSetChart(s->field[f], pStart, pEnd);CHKERRQ(ierr);
633 }
634 PetscFunctionReturn(0);
635 }
636
637 /*@
638 PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
639
640 Not collective
641
642 Input Parameter:
643 . s - the PetscSection
644
645 Output Parameters:
646 . perm - The permutation as an IS
647
648 Level: intermediate
649
650 .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
651 @*/
PetscSectionGetPermutation(PetscSection s,IS * perm)652 PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
653 {
654 PetscFunctionBegin;
655 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
656 if (perm) {PetscValidPointer(perm, 2); *perm = s->perm;}
657 PetscFunctionReturn(0);
658 }
659
660 /*@
661 PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
662
663 Not collective
664
665 Input Parameters:
666 + s - the PetscSection
667 - perm - the permutation of points
668
669 Level: intermediate
670
671 .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
672 @*/
PetscSectionSetPermutation(PetscSection s,IS perm)673 PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
674 {
675 PetscErrorCode ierr;
676
677 PetscFunctionBegin;
678 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
679 if (perm) PetscValidHeaderSpecific(perm, IS_CLASSID, 2);
680 if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
681 if (s->perm != perm) {
682 ierr = ISDestroy(&s->perm);CHKERRQ(ierr);
683 if (perm) {
684 s->perm = perm;
685 ierr = PetscObjectReference((PetscObject) s->perm);CHKERRQ(ierr);
686 }
687 }
688 PetscFunctionReturn(0);
689 }
690
691 /*@
692 PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
693
694 Not collective
695
696 Input Parameter:
697 . s - the PetscSection
698
699 Output Parameter:
700 . pm - the flag for point major ordering
701
702 Level: intermediate
703
704 .seealso: PetscSectionSetPointMajor()
705 @*/
PetscSectionGetPointMajor(PetscSection s,PetscBool * pm)706 PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
707 {
708 PetscFunctionBegin;
709 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
710 PetscValidPointer(pm,2);
711 *pm = s->pointMajor;
712 PetscFunctionReturn(0);
713 }
714
715 /*@
716 PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
717
718 Not collective
719
720 Input Parameters:
721 + s - the PetscSection
722 - pm - the flag for point major ordering
723
724 Not collective
725
726 Level: intermediate
727
728 .seealso: PetscSectionGetPointMajor()
729 @*/
PetscSectionSetPointMajor(PetscSection s,PetscBool pm)730 PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
731 {
732 PetscFunctionBegin;
733 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
734 if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
735 s->pointMajor = pm;
736 PetscFunctionReturn(0);
737 }
738
739 /*@
740 PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
741
742 Not collective
743
744 Input Parameters:
745 + s - the PetscSection
746 - point - the point
747
748 Output Parameter:
749 . numDof - the number of dof
750
751 Level: intermediate
752
753 .seealso: PetscSectionSetDof(), PetscSectionCreate()
754 @*/
PetscSectionGetDof(PetscSection s,PetscInt point,PetscInt * numDof)755 PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
756 {
757 PetscFunctionBeginHot;
758 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
759 PetscValidPointer(numDof, 3);
760 if (PetscDefined(USE_DEBUG)) {
761 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);
762 }
763 *numDof = s->atlasDof[point - s->pStart];
764 PetscFunctionReturn(0);
765 }
766
767 /*@
768 PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
769
770 Not collective
771
772 Input Parameters:
773 + s - the PetscSection
774 . point - the point
775 - numDof - the number of dof
776
777 Level: intermediate
778
779 .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
780 @*/
PetscSectionSetDof(PetscSection s,PetscInt point,PetscInt numDof)781 PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
782 {
783 PetscFunctionBegin;
784 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
785 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);
786 s->atlasDof[point - s->pStart] = numDof;
787 PetscFunctionReturn(0);
788 }
789
790 /*@
791 PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
792
793 Not collective
794
795 Input Parameters:
796 + s - the PetscSection
797 . point - the point
798 - numDof - the number of additional dof
799
800 Level: intermediate
801
802 .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
803 @*/
PetscSectionAddDof(PetscSection s,PetscInt point,PetscInt numDof)804 PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
805 {
806 PetscFunctionBeginHot;
807 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
808 if (PetscDefined(USE_DEBUG)) {
809 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);
810 }
811 s->atlasDof[point - s->pStart] += numDof;
812 PetscFunctionReturn(0);
813 }
814
815 /*@
816 PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
817
818 Not collective
819
820 Input Parameters:
821 + s - the PetscSection
822 . point - the point
823 - field - the field
824
825 Output Parameter:
826 . numDof - the number of dof
827
828 Level: intermediate
829
830 .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
831 @*/
PetscSectionGetFieldDof(PetscSection s,PetscInt point,PetscInt field,PetscInt * numDof)832 PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
833 {
834 PetscErrorCode ierr;
835
836 PetscFunctionBegin;
837 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
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 ierr = PetscSectionGetDof(s->field[field], point, numDof);CHKERRQ(ierr);
840 PetscFunctionReturn(0);
841 }
842
843 /*@
844 PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
845
846 Not collective
847
848 Input Parameters:
849 + s - the PetscSection
850 . point - the point
851 . field - the field
852 - numDof - the number of dof
853
854 Level: intermediate
855
856 .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
857 @*/
PetscSectionSetFieldDof(PetscSection s,PetscInt point,PetscInt field,PetscInt numDof)858 PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
859 {
860 PetscErrorCode ierr;
861
862 PetscFunctionBegin;
863 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
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 ierr = PetscSectionSetDof(s->field[field], point, numDof);CHKERRQ(ierr);
866 PetscFunctionReturn(0);
867 }
868
869 /*@
870 PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
871
872 Not collective
873
874 Input Parameters:
875 + s - the PetscSection
876 . point - the point
877 . field - the field
878 - numDof - the number of dof
879
880 Level: intermediate
881
882 .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
883 @*/
PetscSectionAddFieldDof(PetscSection s,PetscInt point,PetscInt field,PetscInt numDof)884 PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
885 {
886 PetscErrorCode ierr;
887
888 PetscFunctionBegin;
889 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
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 ierr = PetscSectionAddDof(s->field[field], point, numDof);CHKERRQ(ierr);
892 PetscFunctionReturn(0);
893 }
894
895 /*@
896 PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
897
898 Not collective
899
900 Input Parameters:
901 + s - the PetscSection
902 - point - the point
903
904 Output Parameter:
905 . numDof - the number of dof which are fixed by constraints
906
907 Level: intermediate
908
909 .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
910 @*/
PetscSectionGetConstraintDof(PetscSection s,PetscInt point,PetscInt * numDof)911 PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
912 {
913 PetscErrorCode ierr;
914
915 PetscFunctionBegin;
916 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
917 PetscValidPointer(numDof, 3);
918 if (s->bc) {
919 ierr = PetscSectionGetDof(s->bc, point, numDof);CHKERRQ(ierr);
920 } else *numDof = 0;
921 PetscFunctionReturn(0);
922 }
923
924 /*@
925 PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
926
927 Not collective
928
929 Input Parameters:
930 + s - the PetscSection
931 . point - the point
932 - numDof - the number of dof which are fixed by constraints
933
934 Level: intermediate
935
936 .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
937 @*/
PetscSectionSetConstraintDof(PetscSection s,PetscInt point,PetscInt numDof)938 PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
939 {
940 PetscErrorCode ierr;
941
942 PetscFunctionBegin;
943 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
944 if (numDof) {
945 ierr = PetscSectionCheckConstraints_Static(s);CHKERRQ(ierr);
946 ierr = PetscSectionSetDof(s->bc, point, numDof);CHKERRQ(ierr);
947 }
948 PetscFunctionReturn(0);
949 }
950
951 /*@
952 PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
953
954 Not collective
955
956 Input Parameters:
957 + s - the PetscSection
958 . point - the point
959 - numDof - the number of additional dof which are fixed by constraints
960
961 Level: intermediate
962
963 .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
964 @*/
PetscSectionAddConstraintDof(PetscSection s,PetscInt point,PetscInt numDof)965 PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
966 {
967 PetscErrorCode ierr;
968
969 PetscFunctionBegin;
970 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
971 if (numDof) {
972 ierr = PetscSectionCheckConstraints_Static(s);CHKERRQ(ierr);
973 ierr = PetscSectionAddDof(s->bc, point, numDof);CHKERRQ(ierr);
974 }
975 PetscFunctionReturn(0);
976 }
977
978 /*@
979 PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
980
981 Not collective
982
983 Input Parameters:
984 + s - the PetscSection
985 . point - the point
986 - field - the field
987
988 Output Parameter:
989 . numDof - the number of dof which are fixed by constraints
990
991 Level: intermediate
992
993 .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
994 @*/
PetscSectionGetFieldConstraintDof(PetscSection s,PetscInt point,PetscInt field,PetscInt * numDof)995 PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
996 {
997 PetscErrorCode ierr;
998
999 PetscFunctionBegin;
1000 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1001 PetscValidPointer(numDof, 4);
1002 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);
1003 ierr = PetscSectionGetConstraintDof(s->field[field], point, numDof);CHKERRQ(ierr);
1004 PetscFunctionReturn(0);
1005 }
1006
1007 /*@
1008 PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1009
1010 Not collective
1011
1012 Input Parameters:
1013 + s - the PetscSection
1014 . point - the point
1015 . field - the field
1016 - numDof - the number of dof which are fixed by constraints
1017
1018 Level: intermediate
1019
1020 .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1021 @*/
PetscSectionSetFieldConstraintDof(PetscSection s,PetscInt point,PetscInt field,PetscInt numDof)1022 PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1023 {
1024 PetscErrorCode ierr;
1025
1026 PetscFunctionBegin;
1027 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1028 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);
1029 ierr = PetscSectionSetConstraintDof(s->field[field], point, numDof);CHKERRQ(ierr);
1030 PetscFunctionReturn(0);
1031 }
1032
1033 /*@
1034 PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1035
1036 Not collective
1037
1038 Input Parameters:
1039 + s - the PetscSection
1040 . point - the point
1041 . field - the field
1042 - numDof - the number of additional dof which are fixed by constraints
1043
1044 Level: intermediate
1045
1046 .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1047 @*/
PetscSectionAddFieldConstraintDof(PetscSection s,PetscInt point,PetscInt field,PetscInt numDof)1048 PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1049 {
1050 PetscErrorCode ierr;
1051
1052 PetscFunctionBegin;
1053 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1054 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);
1055 ierr = PetscSectionAddConstraintDof(s->field[field], point, numDof);CHKERRQ(ierr);
1056 PetscFunctionReturn(0);
1057 }
1058
1059 /*@
1060 PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1061
1062 Not collective
1063
1064 Input Parameter:
1065 . s - the PetscSection
1066
1067 Level: advanced
1068
1069 .seealso: PetscSectionSetUp(), PetscSectionCreate()
1070 @*/
PetscSectionSetUpBC(PetscSection s)1071 PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1072 {
1073 PetscErrorCode ierr;
1074
1075 PetscFunctionBegin;
1076 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1077 if (s->bc) {
1078 const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
1079
1080 ierr = PetscSectionSetUp(s->bc);CHKERRQ(ierr);
1081 ierr = PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);CHKERRQ(ierr);
1082 }
1083 PetscFunctionReturn(0);
1084 }
1085
1086 /*@
1087 PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1088
1089 Not collective
1090
1091 Input Parameter:
1092 . s - the PetscSection
1093
1094 Level: intermediate
1095
1096 .seealso: PetscSectionCreate()
1097 @*/
PetscSectionSetUp(PetscSection s)1098 PetscErrorCode PetscSectionSetUp(PetscSection s)
1099 {
1100 const PetscInt *pind = NULL;
1101 PetscInt offset = 0, foff, p, f;
1102 PetscErrorCode ierr;
1103
1104 PetscFunctionBegin;
1105 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1106 if (s->setup) PetscFunctionReturn(0);
1107 s->setup = PETSC_TRUE;
1108 /* Set offsets and field offsets for all points */
1109 /* Assume that all fields have the same chart */
1110 if (s->perm) {ierr = ISGetIndices(s->perm, &pind);CHKERRQ(ierr);}
1111 if (s->pointMajor) {
1112 for (p = 0; p < s->pEnd - s->pStart; ++p) {
1113 const PetscInt q = pind ? pind[p] : p;
1114
1115 /* Set point offset */
1116 s->atlasOff[q] = offset;
1117 offset += s->atlasDof[q];
1118 s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1119 /* Set field offset */
1120 for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1121 PetscSection sf = s->field[f];
1122
1123 sf->atlasOff[q] = foff;
1124 foff += sf->atlasDof[q];
1125 }
1126 }
1127 } else {
1128 /* Set field offsets for all points */
1129 for (f = 0; f < s->numFields; ++f) {
1130 PetscSection sf = s->field[f];
1131
1132 for (p = 0; p < s->pEnd - s->pStart; ++p) {
1133 const PetscInt q = pind ? pind[p] : p;
1134
1135 sf->atlasOff[q] = offset;
1136 offset += sf->atlasDof[q];
1137 }
1138 }
1139 /* Disable point offsets since these are unused */
1140 for (p = 0; p < s->pEnd - s->pStart; ++p) {
1141 s->atlasOff[p] = -1;
1142 s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1143 }
1144 }
1145 if (s->perm) {ierr = ISRestoreIndices(s->perm, &pind);CHKERRQ(ierr);}
1146 /* Setup BC sections */
1147 ierr = PetscSectionSetUpBC(s);CHKERRQ(ierr);
1148 for (f = 0; f < s->numFields; ++f) {ierr = PetscSectionSetUpBC(s->field[f]);CHKERRQ(ierr);}
1149 PetscFunctionReturn(0);
1150 }
1151
1152 /*@
1153 PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1154
1155 Not collective
1156
1157 Input Parameters:
1158 . s - the PetscSection
1159
1160 Output Parameter:
1161 . maxDof - the maximum dof
1162
1163 Level: intermediate
1164
1165 .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1166 @*/
PetscSectionGetMaxDof(PetscSection s,PetscInt * maxDof)1167 PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1168 {
1169 PetscFunctionBegin;
1170 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1171 PetscValidPointer(maxDof, 2);
1172 *maxDof = s->maxDof;
1173 PetscFunctionReturn(0);
1174 }
1175
1176 /*@
1177 PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1178
1179 Not collective
1180
1181 Input Parameter:
1182 . s - the PetscSection
1183
1184 Output Parameter:
1185 . size - the size of an array which can hold all the dofs
1186
1187 Level: intermediate
1188
1189 .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1190 @*/
PetscSectionGetStorageSize(PetscSection s,PetscInt * size)1191 PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1192 {
1193 PetscInt p, n = 0;
1194
1195 PetscFunctionBegin;
1196 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1197 PetscValidPointer(size, 2);
1198 for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1199 *size = n;
1200 PetscFunctionReturn(0);
1201 }
1202
1203 /*@
1204 PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1205
1206 Not collective
1207
1208 Input Parameter:
1209 . s - the PetscSection
1210
1211 Output Parameter:
1212 . size - the size of an array which can hold all unconstrained dofs
1213
1214 Level: intermediate
1215
1216 .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1217 @*/
PetscSectionGetConstrainedStorageSize(PetscSection s,PetscInt * size)1218 PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1219 {
1220 PetscInt p, n = 0;
1221
1222 PetscFunctionBegin;
1223 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1224 PetscValidPointer(size, 2);
1225 for (p = 0; p < s->pEnd - s->pStart; ++p) {
1226 const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1227 n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1228 }
1229 *size = n;
1230 PetscFunctionReturn(0);
1231 }
1232
1233 /*@
1234 PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1235 the local section and an SF describing the section point overlap.
1236
1237 Input Parameters:
1238 + s - The PetscSection for the local field layout
1239 . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1240 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1241 - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1242
1243 Output Parameter:
1244 . gsection - The PetscSection for the global field layout
1245
1246 Note: This gives negative sizes and offsets to points not owned by this process
1247
1248 Level: intermediate
1249
1250 .seealso: PetscSectionCreate()
1251 @*/
PetscSectionCreateGlobalSection(PetscSection s,PetscSF sf,PetscBool includeConstraints,PetscBool localOffsets,PetscSection * gsection)1252 PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1253 {
1254 PetscSection gs;
1255 const PetscInt *pind = NULL;
1256 PetscInt *recv = NULL, *neg = NULL;
1257 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1258 PetscErrorCode ierr;
1259
1260 PetscFunctionBegin;
1261 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1262 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2);
1263 PetscValidLogicalCollectiveBool(s, includeConstraints, 3);
1264 PetscValidLogicalCollectiveBool(s, localOffsets, 4);
1265 PetscValidPointer(gsection, 5);
1266 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);CHKERRQ(ierr);
1267 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1268 ierr = PetscSectionSetChart(gs, pStart, pEnd);CHKERRQ(ierr);
1269 ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
1270 nlocal = nroots; /* The local/leaf space matches global/root space */
1271 /* Must allocate for all points visible to SF, which may be more than this section */
1272 if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1273 ierr = PetscSFGetLeafRange(sf, NULL, &maxleaf);CHKERRQ(ierr);
1274 if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1275 if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1276 ierr = PetscMalloc2(nroots,&neg,nlocal,&recv);CHKERRQ(ierr);
1277 ierr = PetscArrayzero(neg,nroots);CHKERRQ(ierr);
1278 }
1279 /* Mark all local points with negative dof */
1280 for (p = pStart; p < pEnd; ++p) {
1281 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1282 ierr = PetscSectionSetDof(gs, p, dof);CHKERRQ(ierr);
1283 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
1284 if (!includeConstraints && cdof > 0) {ierr = PetscSectionSetConstraintDof(gs, p, cdof);CHKERRQ(ierr);}
1285 if (neg) neg[p] = -(dof+1);
1286 }
1287 ierr = PetscSectionSetUpBC(gs);CHKERRQ(ierr);
1288 if (gs->bcIndices) {ierr = PetscArraycpy(gs->bcIndices, s->bcIndices,gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);CHKERRQ(ierr);}
1289 if (nroots >= 0) {
1290 ierr = PetscArrayzero(recv,nlocal);CHKERRQ(ierr);
1291 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, recv);CHKERRQ(ierr);
1292 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, recv);CHKERRQ(ierr);
1293 for (p = pStart; p < pEnd; ++p) {
1294 if (recv[p] < 0) {
1295 gs->atlasDof[p-pStart] = recv[p];
1296 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1297 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);
1298 }
1299 }
1300 }
1301 /* Calculate new sizes, get process offset, and calculate point offsets */
1302 if (s->perm) {ierr = ISGetIndices(s->perm, &pind);CHKERRQ(ierr);}
1303 for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1304 const PetscInt q = pind ? pind[p] : p;
1305
1306 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1307 gs->atlasOff[q] = off;
1308 off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1309 }
1310 if (!localOffsets) {
1311 ierr = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));CHKERRQ(ierr);
1312 globalOff -= off;
1313 }
1314 for (p = pStart, off = 0; p < pEnd; ++p) {
1315 gs->atlasOff[p-pStart] += globalOff;
1316 if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1317 }
1318 if (s->perm) {ierr = ISRestoreIndices(s->perm, &pind);CHKERRQ(ierr);}
1319 /* Put in negative offsets for ghost points */
1320 if (nroots >= 0) {
1321 ierr = PetscArrayzero(recv,nlocal);CHKERRQ(ierr);
1322 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, recv);CHKERRQ(ierr);
1323 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, recv);CHKERRQ(ierr);
1324 for (p = pStart; p < pEnd; ++p) {
1325 if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1326 }
1327 }
1328 ierr = PetscFree2(neg,recv);CHKERRQ(ierr);
1329 ierr = PetscSectionViewFromOptions(gs, NULL, "-global_section_view");CHKERRQ(ierr);
1330 *gsection = gs;
1331 PetscFunctionReturn(0);
1332 }
1333
1334 /*@
1335 PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1336 the local section and an SF describing the section point overlap.
1337
1338 Input Parameters:
1339 + s - The PetscSection for the local field layout
1340 . sf - The SF describing parallel layout of the section points
1341 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1342 . numExcludes - The number of exclusion ranges
1343 - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1344
1345 Output Parameter:
1346 . gsection - The PetscSection for the global field layout
1347
1348 Note: This gives negative sizes and offsets to points not owned by this process
1349
1350 Level: advanced
1351
1352 .seealso: PetscSectionCreate()
1353 @*/
PetscSectionCreateGlobalSectionCensored(PetscSection s,PetscSF sf,PetscBool includeConstraints,PetscInt numExcludes,const PetscInt excludes[],PetscSection * gsection)1354 PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1355 {
1356 const PetscInt *pind = NULL;
1357 PetscInt *neg = NULL, *tmpOff = NULL;
1358 PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1359 PetscErrorCode ierr;
1360
1361 PetscFunctionBegin;
1362 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1363 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2);
1364 PetscValidPointer(gsection, 6);
1365 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);CHKERRQ(ierr);
1366 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1367 ierr = PetscSectionSetChart(*gsection, pStart, pEnd);CHKERRQ(ierr);
1368 ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
1369 if (nroots >= 0) {
1370 if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1371 ierr = PetscCalloc1(nroots, &neg);CHKERRQ(ierr);
1372 if (nroots > pEnd-pStart) {
1373 ierr = PetscCalloc1(nroots, &tmpOff);CHKERRQ(ierr);
1374 } else {
1375 tmpOff = &(*gsection)->atlasDof[-pStart];
1376 }
1377 }
1378 /* Mark ghost points with negative dof */
1379 for (p = pStart; p < pEnd; ++p) {
1380 for (e = 0; e < numExcludes; ++e) {
1381 if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1382 ierr = PetscSectionSetDof(*gsection, p, 0);CHKERRQ(ierr);
1383 break;
1384 }
1385 }
1386 if (e < numExcludes) continue;
1387 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1388 ierr = PetscSectionSetDof(*gsection, p, dof);CHKERRQ(ierr);
1389 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
1390 if (!includeConstraints && cdof > 0) {ierr = PetscSectionSetConstraintDof(*gsection, p, cdof);CHKERRQ(ierr);}
1391 if (neg) neg[p] = -(dof+1);
1392 }
1393 ierr = PetscSectionSetUpBC(*gsection);CHKERRQ(ierr);
1394 if (nroots >= 0) {
1395 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1396 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1397 if (nroots > pEnd - pStart) {
1398 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1399 }
1400 }
1401 /* Calculate new sizes, get proccess offset, and calculate point offsets */
1402 if (s->perm) {ierr = ISGetIndices(s->perm, &pind);CHKERRQ(ierr);}
1403 for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1404 const PetscInt q = pind ? pind[p] : p;
1405
1406 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1407 (*gsection)->atlasOff[q] = off;
1408 off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1409 }
1410 ierr = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));CHKERRQ(ierr);
1411 globalOff -= off;
1412 for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1413 (*gsection)->atlasOff[p] += globalOff;
1414 if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1415 }
1416 if (s->perm) {ierr = ISRestoreIndices(s->perm, &pind);CHKERRQ(ierr);}
1417 /* Put in negative offsets for ghost points */
1418 if (nroots >= 0) {
1419 if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1420 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1421 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr);
1422 if (nroots > pEnd - pStart) {
1423 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1424 }
1425 }
1426 if (nroots >= 0 && nroots > pEnd-pStart) {ierr = PetscFree(tmpOff);CHKERRQ(ierr);}
1427 ierr = PetscFree(neg);CHKERRQ(ierr);
1428 PetscFunctionReturn(0);
1429 }
1430
1431 /*@
1432 PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1433
1434 Collective on comm
1435
1436 Input Parameters:
1437 + comm - The MPI_Comm
1438 - s - The PetscSection
1439
1440 Output Parameter:
1441 . layout - The point layout for the section
1442
1443 Note: This is usually caleld for the default global section.
1444
1445 Level: advanced
1446
1447 .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1448 @*/
PetscSectionGetPointLayout(MPI_Comm comm,PetscSection s,PetscLayout * layout)1449 PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1450 {
1451 PetscInt pStart, pEnd, p, localSize = 0;
1452 PetscErrorCode ierr;
1453
1454 PetscFunctionBegin;
1455 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1456 for (p = pStart; p < pEnd; ++p) {
1457 PetscInt dof;
1458
1459 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1460 if (dof > 0) ++localSize;
1461 }
1462 ierr = PetscLayoutCreate(comm, layout);CHKERRQ(ierr);
1463 ierr = PetscLayoutSetLocalSize(*layout, localSize);CHKERRQ(ierr);
1464 ierr = PetscLayoutSetBlockSize(*layout, 1);CHKERRQ(ierr);
1465 ierr = PetscLayoutSetUp(*layout);CHKERRQ(ierr);
1466 PetscFunctionReturn(0);
1467 }
1468
1469 /*@
1470 PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1471
1472 Collective on comm
1473
1474 Input Parameters:
1475 + comm - The MPI_Comm
1476 - s - The PetscSection
1477
1478 Output Parameter:
1479 . layout - The dof layout for the section
1480
1481 Note: This is usually called for the default global section.
1482
1483 Level: advanced
1484
1485 .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1486 @*/
PetscSectionGetValueLayout(MPI_Comm comm,PetscSection s,PetscLayout * layout)1487 PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1488 {
1489 PetscInt pStart, pEnd, p, localSize = 0;
1490 PetscErrorCode ierr;
1491
1492 PetscFunctionBegin;
1493 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2);
1494 PetscValidPointer(layout, 3);
1495 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1496 for (p = pStart; p < pEnd; ++p) {
1497 PetscInt dof,cdof;
1498
1499 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1500 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
1501 if (dof-cdof > 0) localSize += dof-cdof;
1502 }
1503 ierr = PetscLayoutCreate(comm, layout);CHKERRQ(ierr);
1504 ierr = PetscLayoutSetLocalSize(*layout, localSize);CHKERRQ(ierr);
1505 ierr = PetscLayoutSetBlockSize(*layout, 1);CHKERRQ(ierr);
1506 ierr = PetscLayoutSetUp(*layout);CHKERRQ(ierr);
1507 PetscFunctionReturn(0);
1508 }
1509
1510 /*@
1511 PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1512
1513 Not collective
1514
1515 Input Parameters:
1516 + s - the PetscSection
1517 - point - the point
1518
1519 Output Parameter:
1520 . offset - the offset
1521
1522 Level: intermediate
1523
1524 .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1525 @*/
PetscSectionGetOffset(PetscSection s,PetscInt point,PetscInt * offset)1526 PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1527 {
1528 PetscFunctionBegin;
1529 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1530 PetscValidPointer(offset, 3);
1531 if (PetscDefined(USE_DEBUG)) {
1532 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);
1533 }
1534 *offset = s->atlasOff[point - s->pStart];
1535 PetscFunctionReturn(0);
1536 }
1537
1538 /*@
1539 PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1540
1541 Not collective
1542
1543 Input Parameters:
1544 + s - the PetscSection
1545 . point - the point
1546 - offset - the offset
1547
1548 Note: The user usually does not call this function, but uses PetscSectionSetUp()
1549
1550 Level: intermediate
1551
1552 .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1553 @*/
PetscSectionSetOffset(PetscSection s,PetscInt point,PetscInt offset)1554 PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1555 {
1556 PetscFunctionBegin;
1557 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1558 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);
1559 s->atlasOff[point - s->pStart] = offset;
1560 PetscFunctionReturn(0);
1561 }
1562
1563 /*@
1564 PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1565
1566 Not collective
1567
1568 Input Parameters:
1569 + s - the PetscSection
1570 . point - the point
1571 - field - the field
1572
1573 Output Parameter:
1574 . offset - the offset
1575
1576 Level: intermediate
1577
1578 .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1579 @*/
PetscSectionGetFieldOffset(PetscSection s,PetscInt point,PetscInt field,PetscInt * offset)1580 PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1581 {
1582 PetscErrorCode ierr;
1583
1584 PetscFunctionBegin;
1585 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1586 PetscValidPointer(offset, 4);
1587 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);
1588 ierr = PetscSectionGetOffset(s->field[field], point, offset);CHKERRQ(ierr);
1589 PetscFunctionReturn(0);
1590 }
1591
1592 /*@
1593 PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1594
1595 Not collective
1596
1597 Input Parameters:
1598 + s - the PetscSection
1599 . point - the point
1600 . field - the field
1601 - offset - the offset
1602
1603 Note: The user usually does not call this function, but uses PetscSectionSetUp()
1604
1605 Level: intermediate
1606
1607 .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1608 @*/
PetscSectionSetFieldOffset(PetscSection s,PetscInt point,PetscInt field,PetscInt offset)1609 PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1610 {
1611 PetscErrorCode ierr;
1612
1613 PetscFunctionBegin;
1614 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1615 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);
1616 ierr = PetscSectionSetOffset(s->field[field], point, offset);CHKERRQ(ierr);
1617 PetscFunctionReturn(0);
1618 }
1619
1620 /*@
1621 PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1622
1623 Not collective
1624
1625 Input Parameters:
1626 + s - the PetscSection
1627 . point - the point
1628 - field - the field
1629
1630 Output Parameter:
1631 . offset - the offset
1632
1633 Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1634 this point, what number is the first dof with this field.
1635
1636 Level: advanced
1637
1638 .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1639 @*/
PetscSectionGetFieldPointOffset(PetscSection s,PetscInt point,PetscInt field,PetscInt * offset)1640 PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1641 {
1642 PetscInt off, foff;
1643 PetscErrorCode ierr;
1644
1645 PetscFunctionBegin;
1646 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1647 PetscValidPointer(offset, 4);
1648 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);
1649 ierr = PetscSectionGetOffset(s, point, &off);CHKERRQ(ierr);
1650 ierr = PetscSectionGetOffset(s->field[field], point, &foff);CHKERRQ(ierr);
1651 *offset = foff - off;
1652 PetscFunctionReturn(0);
1653 }
1654
1655 /*@
1656 PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1657
1658 Not collective
1659
1660 Input Parameter:
1661 . s - the PetscSection
1662
1663 Output Parameters:
1664 + start - the minimum offset
1665 - end - one more than the maximum offset
1666
1667 Level: intermediate
1668
1669 .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1670 @*/
PetscSectionGetOffsetRange(PetscSection s,PetscInt * start,PetscInt * end)1671 PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1672 {
1673 PetscInt os = 0, oe = 0, pStart, pEnd, p;
1674 PetscErrorCode ierr;
1675
1676 PetscFunctionBegin;
1677 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1678 if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1679 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1680 for (p = 0; p < pEnd-pStart; ++p) {
1681 PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1682
1683 if (off >= 0) {
1684 os = PetscMin(os, off);
1685 oe = PetscMax(oe, off+dof);
1686 }
1687 }
1688 if (start) *start = os;
1689 if (end) *end = oe;
1690 PetscFunctionReturn(0);
1691 }
1692
1693 /*@
1694 PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1695
1696 Collective on s
1697
1698 Input Parameter:
1699 + s - the PetscSection
1700 . len - the number of subfields
1701 - fields - the subfield numbers
1702
1703 Output Parameter:
1704 . subs - the subsection
1705
1706 Note: The section offsets now refer to a new, smaller vector.
1707
1708 Level: advanced
1709
1710 .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1711 @*/
PetscSectionCreateSubsection(PetscSection s,PetscInt len,const PetscInt fields[],PetscSection * subs)1712 PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1713 {
1714 PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1715 PetscErrorCode ierr;
1716
1717 PetscFunctionBegin;
1718 if (!len) PetscFunctionReturn(0);
1719 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1720 PetscValidPointer(fields, 3);
1721 PetscValidPointer(subs, 4);
1722 ierr = PetscSectionGetNumFields(s, &nF);CHKERRQ(ierr);
1723 if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1724 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);CHKERRQ(ierr);
1725 ierr = PetscSectionSetNumFields(*subs, len);CHKERRQ(ierr);
1726 for (f = 0; f < len; ++f) {
1727 const char *name = NULL;
1728 PetscInt numComp = 0;
1729
1730 ierr = PetscSectionGetFieldName(s, fields[f], &name);CHKERRQ(ierr);
1731 ierr = PetscSectionSetFieldName(*subs, f, name);CHKERRQ(ierr);
1732 ierr = PetscSectionGetFieldComponents(s, fields[f], &numComp);CHKERRQ(ierr);
1733 ierr = PetscSectionSetFieldComponents(*subs, f, numComp);CHKERRQ(ierr);
1734 for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1735 ierr = PetscSectionGetComponentName(s, fields[f], c, &name);CHKERRQ(ierr);
1736 ierr = PetscSectionSetComponentName(*subs, f, c, name);CHKERRQ(ierr);
1737 }
1738 }
1739 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1740 ierr = PetscSectionSetChart(*subs, pStart, pEnd);CHKERRQ(ierr);
1741 for (p = pStart; p < pEnd; ++p) {
1742 PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1743
1744 for (f = 0; f < len; ++f) {
1745 ierr = PetscSectionGetFieldDof(s, p, fields[f], &fdof);CHKERRQ(ierr);
1746 ierr = PetscSectionSetFieldDof(*subs, p, f, fdof);CHKERRQ(ierr);
1747 ierr = PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);CHKERRQ(ierr);
1748 if (cfdof) {ierr = PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);CHKERRQ(ierr);}
1749 dof += fdof;
1750 cdof += cfdof;
1751 }
1752 ierr = PetscSectionSetDof(*subs, p, dof);CHKERRQ(ierr);
1753 if (cdof) {ierr = PetscSectionSetConstraintDof(*subs, p, cdof);CHKERRQ(ierr);}
1754 maxCdof = PetscMax(cdof, maxCdof);
1755 }
1756 ierr = PetscSectionSetUp(*subs);CHKERRQ(ierr);
1757 if (maxCdof) {
1758 PetscInt *indices;
1759
1760 ierr = PetscMalloc1(maxCdof, &indices);CHKERRQ(ierr);
1761 for (p = pStart; p < pEnd; ++p) {
1762 PetscInt cdof;
1763
1764 ierr = PetscSectionGetConstraintDof(*subs, p, &cdof);CHKERRQ(ierr);
1765 if (cdof) {
1766 const PetscInt *oldIndices = NULL;
1767 PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1768
1769 for (f = 0; f < len; ++f) {
1770 PetscInt oldFoff = 0, oldf;
1771
1772 ierr = PetscSectionGetFieldDof(s, p, fields[f], &fdof);CHKERRQ(ierr);
1773 ierr = PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);CHKERRQ(ierr);
1774 ierr = PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);CHKERRQ(ierr);
1775 /* This can be sped up if we assume sorted fields */
1776 for (oldf = 0; oldf < fields[f]; ++oldf) {
1777 PetscInt oldfdof = 0;
1778 ierr = PetscSectionGetFieldDof(s, p, oldf, &oldfdof);CHKERRQ(ierr);
1779 oldFoff += oldfdof;
1780 }
1781 for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1782 ierr = PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);CHKERRQ(ierr);
1783 numConst += cfdof;
1784 fOff += fdof;
1785 }
1786 ierr = PetscSectionSetConstraintIndices(*subs, p, indices);CHKERRQ(ierr);
1787 }
1788 }
1789 ierr = PetscFree(indices);CHKERRQ(ierr);
1790 }
1791 PetscFunctionReturn(0);
1792 }
1793
1794 /*@
1795 PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1796
1797 Collective on s
1798
1799 Input Parameters:
1800 + s - the input sections
1801 - len - the number of input sections
1802
1803 Output Parameter:
1804 . supers - the supersection
1805
1806 Note: The section offsets now refer to a new, larger vector.
1807
1808 Level: advanced
1809
1810 .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1811 @*/
PetscSectionCreateSupersection(PetscSection s[],PetscInt len,PetscSection * supers)1812 PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1813 {
1814 PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1815 PetscErrorCode ierr;
1816
1817 PetscFunctionBegin;
1818 if (!len) PetscFunctionReturn(0);
1819 for (i = 0; i < len; ++i) {
1820 PetscInt nf, pStarti, pEndi;
1821
1822 ierr = PetscSectionGetNumFields(s[i], &nf);CHKERRQ(ierr);
1823 ierr = PetscSectionGetChart(s[i], &pStarti, &pEndi);CHKERRQ(ierr);
1824 pStart = PetscMin(pStart, pStarti);
1825 pEnd = PetscMax(pEnd, pEndi);
1826 Nf += nf;
1827 }
1828 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);CHKERRQ(ierr);
1829 ierr = PetscSectionSetNumFields(*supers, Nf);CHKERRQ(ierr);
1830 for (i = 0, f = 0; i < len; ++i) {
1831 PetscInt nf, fi, ci;
1832
1833 ierr = PetscSectionGetNumFields(s[i], &nf);CHKERRQ(ierr);
1834 for (fi = 0; fi < nf; ++fi, ++f) {
1835 const char *name = NULL;
1836 PetscInt numComp = 0;
1837
1838 ierr = PetscSectionGetFieldName(s[i], fi, &name);CHKERRQ(ierr);
1839 ierr = PetscSectionSetFieldName(*supers, f, name);CHKERRQ(ierr);
1840 ierr = PetscSectionGetFieldComponents(s[i], fi, &numComp);CHKERRQ(ierr);
1841 ierr = PetscSectionSetFieldComponents(*supers, f, numComp);CHKERRQ(ierr);
1842 for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1843 ierr = PetscSectionGetComponentName(s[i], fi, ci, &name);CHKERRQ(ierr);
1844 ierr = PetscSectionSetComponentName(*supers, f, ci, name);CHKERRQ(ierr);
1845 }
1846 }
1847 }
1848 ierr = PetscSectionSetChart(*supers, pStart, pEnd);CHKERRQ(ierr);
1849 for (p = pStart; p < pEnd; ++p) {
1850 PetscInt dof = 0, cdof = 0;
1851
1852 for (i = 0, f = 0; i < len; ++i) {
1853 PetscInt nf, fi, pStarti, pEndi;
1854 PetscInt fdof = 0, cfdof = 0;
1855
1856 ierr = PetscSectionGetNumFields(s[i], &nf);CHKERRQ(ierr);
1857 ierr = PetscSectionGetChart(s[i], &pStarti, &pEndi);CHKERRQ(ierr);
1858 if ((p < pStarti) || (p >= pEndi)) continue;
1859 for (fi = 0; fi < nf; ++fi, ++f) {
1860 ierr = PetscSectionGetFieldDof(s[i], p, fi, &fdof);CHKERRQ(ierr);
1861 ierr = PetscSectionAddFieldDof(*supers, p, f, fdof);CHKERRQ(ierr);
1862 ierr = PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);CHKERRQ(ierr);
1863 if (cfdof) {ierr = PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);CHKERRQ(ierr);}
1864 dof += fdof;
1865 cdof += cfdof;
1866 }
1867 }
1868 ierr = PetscSectionSetDof(*supers, p, dof);CHKERRQ(ierr);
1869 if (cdof) {ierr = PetscSectionSetConstraintDof(*supers, p, cdof);CHKERRQ(ierr);}
1870 maxCdof = PetscMax(cdof, maxCdof);
1871 }
1872 ierr = PetscSectionSetUp(*supers);CHKERRQ(ierr);
1873 if (maxCdof) {
1874 PetscInt *indices;
1875
1876 ierr = PetscMalloc1(maxCdof, &indices);CHKERRQ(ierr);
1877 for (p = pStart; p < pEnd; ++p) {
1878 PetscInt cdof;
1879
1880 ierr = PetscSectionGetConstraintDof(*supers, p, &cdof);CHKERRQ(ierr);
1881 if (cdof) {
1882 PetscInt dof, numConst = 0, fOff = 0;
1883
1884 for (i = 0, f = 0; i < len; ++i) {
1885 const PetscInt *oldIndices = NULL;
1886 PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1887
1888 ierr = PetscSectionGetNumFields(s[i], &nf);CHKERRQ(ierr);
1889 ierr = PetscSectionGetChart(s[i], &pStarti, &pEndi);CHKERRQ(ierr);
1890 if ((p < pStarti) || (p >= pEndi)) continue;
1891 for (fi = 0; fi < nf; ++fi, ++f) {
1892 ierr = PetscSectionGetFieldDof(s[i], p, fi, &fdof);CHKERRQ(ierr);
1893 ierr = PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);CHKERRQ(ierr);
1894 ierr = PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);CHKERRQ(ierr);
1895 for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1896 ierr = PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);CHKERRQ(ierr);
1897 numConst += cfdof;
1898 }
1899 ierr = PetscSectionGetDof(s[i], p, &dof);CHKERRQ(ierr);
1900 fOff += dof;
1901 }
1902 ierr = PetscSectionSetConstraintIndices(*supers, p, indices);CHKERRQ(ierr);
1903 }
1904 }
1905 ierr = PetscFree(indices);CHKERRQ(ierr);
1906 }
1907 PetscFunctionReturn(0);
1908 }
1909
1910 /*@
1911 PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1912
1913 Collective on s
1914
1915 Input Parameters:
1916 + s - the PetscSection
1917 - subpointMap - a sorted list of points in the original mesh which are in the submesh
1918
1919 Output Parameter:
1920 . subs - the subsection
1921
1922 Note: The section offsets now refer to a new, smaller vector.
1923
1924 Level: advanced
1925
1926 .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1927 @*/
PetscSectionCreateSubmeshSection(PetscSection s,IS subpointMap,PetscSection * subs)1928 PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1929 {
1930 const PetscInt *points = NULL, *indices = NULL;
1931 PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;
1932 PetscErrorCode ierr;
1933
1934 PetscFunctionBegin;
1935 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
1936 PetscValidHeaderSpecific(subpointMap, IS_CLASSID, 2);
1937 PetscValidPointer(subs, 3);
1938 ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr);
1939 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);CHKERRQ(ierr);
1940 if (numFields) {ierr = PetscSectionSetNumFields(*subs, numFields);CHKERRQ(ierr);}
1941 for (f = 0; f < numFields; ++f) {
1942 const char *name = NULL;
1943 PetscInt numComp = 0;
1944
1945 ierr = PetscSectionGetFieldName(s, f, &name);CHKERRQ(ierr);
1946 ierr = PetscSectionSetFieldName(*subs, f, name);CHKERRQ(ierr);
1947 ierr = PetscSectionGetFieldComponents(s, f, &numComp);CHKERRQ(ierr);
1948 ierr = PetscSectionSetFieldComponents(*subs, f, numComp);CHKERRQ(ierr);
1949 for (c = 0; c < s->numFieldComponents[f]; ++c) {
1950 ierr = PetscSectionGetComponentName(s, f, c, &name);CHKERRQ(ierr);
1951 ierr = PetscSectionSetComponentName(*subs, f, c, name);CHKERRQ(ierr);
1952 }
1953 }
1954 /* For right now, we do not try to squeeze the subchart */
1955 if (subpointMap) {
1956 ierr = ISGetSize(subpointMap, &numSubpoints);CHKERRQ(ierr);
1957 ierr = ISGetIndices(subpointMap, &points);CHKERRQ(ierr);
1958 }
1959 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
1960 ierr = PetscSectionSetChart(*subs, 0, numSubpoints);CHKERRQ(ierr);
1961 for (p = pStart; p < pEnd; ++p) {
1962 PetscInt dof, cdof, fdof = 0, cfdof = 0;
1963
1964 ierr = PetscFindInt(p, numSubpoints, points, &subp);CHKERRQ(ierr);
1965 if (subp < 0) continue;
1966 for (f = 0; f < numFields; ++f) {
1967 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr);
1968 ierr = PetscSectionSetFieldDof(*subs, subp, f, fdof);CHKERRQ(ierr);
1969 ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);CHKERRQ(ierr);
1970 if (cfdof) {ierr = PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);CHKERRQ(ierr);}
1971 }
1972 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
1973 ierr = PetscSectionSetDof(*subs, subp, dof);CHKERRQ(ierr);
1974 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
1975 if (cdof) {ierr = PetscSectionSetConstraintDof(*subs, subp, cdof);CHKERRQ(ierr);}
1976 }
1977 ierr = PetscSectionSetUp(*subs);CHKERRQ(ierr);
1978 /* Change offsets to original offsets */
1979 for (p = pStart; p < pEnd; ++p) {
1980 PetscInt off, foff = 0;
1981
1982 ierr = PetscFindInt(p, numSubpoints, points, &subp);CHKERRQ(ierr);
1983 if (subp < 0) continue;
1984 for (f = 0; f < numFields; ++f) {
1985 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr);
1986 ierr = PetscSectionSetFieldOffset(*subs, subp, f, foff);CHKERRQ(ierr);
1987 }
1988 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr);
1989 ierr = PetscSectionSetOffset(*subs, subp, off);CHKERRQ(ierr);
1990 }
1991 /* Copy constraint indices */
1992 for (subp = 0; subp < numSubpoints; ++subp) {
1993 PetscInt cdof;
1994
1995 ierr = PetscSectionGetConstraintDof(*subs, subp, &cdof);CHKERRQ(ierr);
1996 if (cdof) {
1997 for (f = 0; f < numFields; ++f) {
1998 ierr = PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);CHKERRQ(ierr);
1999 ierr = PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);CHKERRQ(ierr);
2000 }
2001 ierr = PetscSectionGetConstraintIndices(s, points[subp], &indices);CHKERRQ(ierr);
2002 ierr = PetscSectionSetConstraintIndices(*subs, subp, indices);CHKERRQ(ierr);
2003 }
2004 }
2005 if (subpointMap) {ierr = ISRestoreIndices(subpointMap, &points);CHKERRQ(ierr);}
2006 PetscFunctionReturn(0);
2007 }
2008
PetscSectionView_ASCII(PetscSection s,PetscViewer viewer)2009 static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2010 {
2011 PetscInt p;
2012 PetscMPIInt rank;
2013 PetscErrorCode ierr;
2014
2015 PetscFunctionBegin;
2016 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);CHKERRQ(ierr);
2017 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
2018 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);CHKERRQ(ierr);
2019 for (p = 0; p < s->pEnd - s->pStart; ++p) {
2020 if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2021 PetscInt b;
2022
2023 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);CHKERRQ(ierr);
2024 if (s->bcIndices) {
2025 for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2026 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);CHKERRQ(ierr);
2027 }
2028 }
2029 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "\n");CHKERRQ(ierr);
2030 } else {
2031 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);CHKERRQ(ierr);
2032 }
2033 }
2034 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
2035 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
2036 if (s->sym) {
2037 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
2038 ierr = PetscSectionSymView(s->sym,viewer);CHKERRQ(ierr);
2039 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
2040 }
2041 PetscFunctionReturn(0);
2042 }
2043
2044 /*@C
2045 PetscSectionViewFromOptions - View from Options
2046
2047 Collective on PetscSection
2048
2049 Input Parameters:
2050 + A - the PetscSection object to view
2051 . obj - Optional object
2052 - name - command line option
2053
2054 Level: intermediate
2055 .seealso: PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2056 @*/
PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])2057 PetscErrorCode PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2058 {
2059 PetscErrorCode ierr;
2060
2061 PetscFunctionBegin;
2062 PetscValidHeaderSpecific(A,PETSC_SECTION_CLASSID,1);
2063 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr);
2064 PetscFunctionReturn(0);
2065 }
2066
2067 /*@C
2068 PetscSectionView - Views a PetscSection
2069
2070 Collective on PetscSection
2071
2072 Input Parameters:
2073 + s - the PetscSection object to view
2074 - v - the viewer
2075
2076 Level: beginner
2077
2078 .seealso PetscSectionCreate(), PetscSectionDestroy()
2079 @*/
PetscSectionView(PetscSection s,PetscViewer viewer)2080 PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2081 {
2082 PetscBool isascii;
2083 PetscInt f;
2084 PetscErrorCode ierr;
2085
2086 PetscFunctionBegin;
2087 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2088 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);CHKERRQ(ierr);}
2089 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
2090 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
2091 if (isascii) {
2092 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);CHKERRQ(ierr);
2093 if (s->numFields) {
2094 ierr = PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);CHKERRQ(ierr);
2095 for (f = 0; f < s->numFields; ++f) {
2096 ierr = PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);CHKERRQ(ierr);
2097 ierr = PetscSectionView_ASCII(s->field[f], viewer);CHKERRQ(ierr);
2098 }
2099 } else {
2100 ierr = PetscSectionView_ASCII(s, viewer);CHKERRQ(ierr);
2101 }
2102 }
2103 PetscFunctionReturn(0);
2104 }
2105
PetscSectionResetClosurePermutation(PetscSection section)2106 static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2107 {
2108 PetscErrorCode ierr;
2109 PetscSectionClosurePermVal clVal;
2110
2111 PetscFunctionBegin;
2112 if (!section->clHash) PetscFunctionReturn(0);
2113 kh_foreach_value(section->clHash, clVal, {
2114 ierr = PetscFree(clVal.perm);CHKERRQ(ierr);
2115 ierr = PetscFree(clVal.invPerm);CHKERRQ(ierr);
2116 });
2117 kh_destroy(ClPerm, section->clHash);
2118 section->clHash = NULL;
2119 PetscFunctionReturn(0);
2120 }
2121
2122 /*@
2123 PetscSectionReset - Frees all section data.
2124
2125 Not collective
2126
2127 Input Parameters:
2128 . s - the PetscSection
2129
2130 Level: beginner
2131
2132 .seealso: PetscSection, PetscSectionCreate()
2133 @*/
PetscSectionReset(PetscSection s)2134 PetscErrorCode PetscSectionReset(PetscSection s)
2135 {
2136 PetscInt f, c;
2137 PetscErrorCode ierr;
2138
2139 PetscFunctionBegin;
2140 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2141 for (f = 0; f < s->numFields; ++f) {
2142 ierr = PetscSectionDestroy(&s->field[f]);CHKERRQ(ierr);
2143 ierr = PetscFree(s->fieldNames[f]);CHKERRQ(ierr);
2144 for (c = 0; c < s->numFieldComponents[f]; ++c)
2145 ierr = PetscFree(s->compNames[f][c]);CHKERRQ(ierr);
2146 ierr = PetscFree(s->compNames[f]);CHKERRQ(ierr);
2147 }
2148 ierr = PetscFree(s->numFieldComponents);CHKERRQ(ierr);
2149 ierr = PetscFree(s->fieldNames);CHKERRQ(ierr);
2150 ierr = PetscFree(s->compNames);CHKERRQ(ierr);
2151 ierr = PetscFree(s->field);CHKERRQ(ierr);
2152 ierr = PetscSectionDestroy(&s->bc);CHKERRQ(ierr);
2153 ierr = PetscFree(s->bcIndices);CHKERRQ(ierr);
2154 ierr = PetscFree2(s->atlasDof, s->atlasOff);CHKERRQ(ierr);
2155 ierr = PetscSectionDestroy(&s->clSection);CHKERRQ(ierr);
2156 ierr = ISDestroy(&s->clPoints);CHKERRQ(ierr);
2157 ierr = ISDestroy(&s->perm);CHKERRQ(ierr);
2158 ierr = PetscSectionResetClosurePermutation(s);CHKERRQ(ierr);
2159 ierr = PetscSectionSymDestroy(&s->sym);CHKERRQ(ierr);
2160 ierr = PetscSectionDestroy(&s->clSection);CHKERRQ(ierr);
2161 ierr = ISDestroy(&s->clPoints);CHKERRQ(ierr);
2162
2163 s->pStart = -1;
2164 s->pEnd = -1;
2165 s->maxDof = 0;
2166 s->setup = PETSC_FALSE;
2167 s->numFields = 0;
2168 s->clObj = NULL;
2169 PetscFunctionReturn(0);
2170 }
2171
2172 /*@
2173 PetscSectionDestroy - Frees a section object and frees its range if that exists.
2174
2175 Not collective
2176
2177 Input Parameters:
2178 . s - the PetscSection
2179
2180 Level: beginner
2181
2182 .seealso: PetscSection, PetscSectionCreate()
2183 @*/
PetscSectionDestroy(PetscSection * s)2184 PetscErrorCode PetscSectionDestroy(PetscSection *s)
2185 {
2186 PetscErrorCode ierr;
2187
2188 PetscFunctionBegin;
2189 if (!*s) PetscFunctionReturn(0);
2190 PetscValidHeaderSpecific(*s,PETSC_SECTION_CLASSID,1);
2191 if (--((PetscObject)(*s))->refct > 0) {
2192 *s = NULL;
2193 PetscFunctionReturn(0);
2194 }
2195 ierr = PetscSectionReset(*s);CHKERRQ(ierr);
2196 ierr = PetscHeaderDestroy(s);CHKERRQ(ierr);
2197 PetscFunctionReturn(0);
2198 }
2199
VecIntGetValuesSection(PetscInt * baseArray,PetscSection s,PetscInt point,const PetscInt ** values)2200 PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2201 {
2202 const PetscInt p = point - s->pStart;
2203
2204 PetscFunctionBegin;
2205 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2);
2206 *values = &baseArray[s->atlasOff[p]];
2207 PetscFunctionReturn(0);
2208 }
2209
VecIntSetValuesSection(PetscInt * baseArray,PetscSection s,PetscInt point,const PetscInt values[],InsertMode mode)2210 PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2211 {
2212 PetscInt *array;
2213 const PetscInt p = point - s->pStart;
2214 const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2215 PetscInt cDim = 0;
2216 PetscErrorCode ierr;
2217
2218 PetscFunctionBegin;
2219 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2);
2220 ierr = PetscSectionGetConstraintDof(s, p, &cDim);CHKERRQ(ierr);
2221 array = &baseArray[s->atlasOff[p]];
2222 if (!cDim) {
2223 if (orientation >= 0) {
2224 const PetscInt dim = s->atlasDof[p];
2225 PetscInt i;
2226
2227 if (mode == INSERT_VALUES) {
2228 for (i = 0; i < dim; ++i) array[i] = values[i];
2229 } else {
2230 for (i = 0; i < dim; ++i) array[i] += values[i];
2231 }
2232 } else {
2233 PetscInt offset = 0;
2234 PetscInt j = -1, field, i;
2235
2236 for (field = 0; field < s->numFields; ++field) {
2237 const PetscInt dim = s->field[field]->atlasDof[p];
2238
2239 for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2240 offset += dim;
2241 }
2242 }
2243 } else {
2244 if (orientation >= 0) {
2245 const PetscInt dim = s->atlasDof[p];
2246 PetscInt cInd = 0, i;
2247 const PetscInt *cDof;
2248
2249 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr);
2250 if (mode == INSERT_VALUES) {
2251 for (i = 0; i < dim; ++i) {
2252 if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2253 array[i] = values[i];
2254 }
2255 } else {
2256 for (i = 0; i < dim; ++i) {
2257 if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2258 array[i] += values[i];
2259 }
2260 }
2261 } else {
2262 const PetscInt *cDof;
2263 PetscInt offset = 0;
2264 PetscInt cOffset = 0;
2265 PetscInt j = 0, field;
2266
2267 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr);
2268 for (field = 0; field < s->numFields; ++field) {
2269 const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2270 const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2271 const PetscInt sDim = dim - tDim;
2272 PetscInt cInd = 0, i,k;
2273
2274 for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2275 if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2276 array[j] = values[k];
2277 }
2278 offset += dim;
2279 cOffset += dim - tDim;
2280 }
2281 }
2282 }
2283 PetscFunctionReturn(0);
2284 }
2285
2286 /*@C
2287 PetscSectionHasConstraints - Determine whether a section has constrained dofs
2288
2289 Not collective
2290
2291 Input Parameter:
2292 . s - The PetscSection
2293
2294 Output Parameter:
2295 . hasConstraints - flag indicating that the section has constrained dofs
2296
2297 Level: intermediate
2298
2299 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2300 @*/
PetscSectionHasConstraints(PetscSection s,PetscBool * hasConstraints)2301 PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2302 {
2303 PetscFunctionBegin;
2304 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2305 PetscValidPointer(hasConstraints, 2);
2306 *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2307 PetscFunctionReturn(0);
2308 }
2309
2310 /*@C
2311 PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2312
2313 Not collective
2314
2315 Input Parameters:
2316 + s - The PetscSection
2317 - point - The point
2318
2319 Output Parameter:
2320 . indices - The constrained dofs
2321
2322 Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2323
2324 Level: intermediate
2325
2326 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2327 @*/
PetscSectionGetConstraintIndices(PetscSection s,PetscInt point,const PetscInt ** indices)2328 PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2329 {
2330 PetscErrorCode ierr;
2331
2332 PetscFunctionBegin;
2333 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2334 if (s->bc) {
2335 ierr = VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);CHKERRQ(ierr);
2336 } else *indices = NULL;
2337 PetscFunctionReturn(0);
2338 }
2339
2340 /*@C
2341 PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2342
2343 Not collective
2344
2345 Input Parameters:
2346 + s - The PetscSection
2347 . point - The point
2348 - indices - The constrained dofs
2349
2350 Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2351
2352 Level: intermediate
2353
2354 .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2355 @*/
PetscSectionSetConstraintIndices(PetscSection s,PetscInt point,const PetscInt indices[])2356 PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2357 {
2358 PetscErrorCode ierr;
2359
2360 PetscFunctionBegin;
2361 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2362 if (s->bc) {
2363 ierr = VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);CHKERRQ(ierr);
2364 }
2365 PetscFunctionReturn(0);
2366 }
2367
2368 /*@C
2369 PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2370
2371 Not collective
2372
2373 Input Parameters:
2374 + s - The PetscSection
2375 . field - The field number
2376 - point - The point
2377
2378 Output Parameter:
2379 . indices - The constrained dofs sorted in ascending order
2380
2381 Notes:
2382 The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().
2383
2384 Fortran Note:
2385 In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2386
2387 Level: intermediate
2388
2389 .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2390 @*/
PetscSectionGetFieldConstraintIndices(PetscSection s,PetscInt point,PetscInt field,const PetscInt ** indices)2391 PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2392 {
2393 PetscErrorCode ierr;
2394
2395 PetscFunctionBegin;
2396 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2397 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);
2398 ierr = PetscSectionGetConstraintIndices(s->field[field], point, indices);CHKERRQ(ierr);
2399 PetscFunctionReturn(0);
2400 }
2401
2402 /*@C
2403 PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2404
2405 Not collective
2406
2407 Input Parameters:
2408 + s - The PetscSection
2409 . point - The point
2410 . field - The field number
2411 - indices - The constrained dofs
2412
2413 Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2414
2415 Level: intermediate
2416
2417 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2418 @*/
PetscSectionSetFieldConstraintIndices(PetscSection s,PetscInt point,PetscInt field,const PetscInt indices[])2419 PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2420 {
2421 PetscErrorCode ierr;
2422
2423 PetscFunctionBegin;
2424 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
2425 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);
2426 ierr = PetscSectionSetConstraintIndices(s->field[field], point, indices);CHKERRQ(ierr);
2427 PetscFunctionReturn(0);
2428 }
2429
2430 /*@
2431 PetscSectionPermute - Reorder the section according to the input point permutation
2432
2433 Collective on PetscSection
2434
2435 Input Parameter:
2436 + section - The PetscSection object
2437 - perm - The point permutation, old point p becomes new point perm[p]
2438
2439 Output Parameter:
2440 . sectionNew - The permuted PetscSection
2441
2442 Level: intermediate
2443
2444 .seealso: MatPermute()
2445 @*/
PetscSectionPermute(PetscSection section,IS permutation,PetscSection * sectionNew)2446 PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2447 {
2448 PetscSection s = section, sNew;
2449 const PetscInt *perm;
2450 PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2451 PetscErrorCode ierr;
2452
2453 PetscFunctionBegin;
2454 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1);
2455 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2);
2456 PetscValidPointer(sectionNew, 3);
2457 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);CHKERRQ(ierr);
2458 ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr);
2459 if (numFields) {ierr = PetscSectionSetNumFields(sNew, numFields);CHKERRQ(ierr);}
2460 for (f = 0; f < numFields; ++f) {
2461 const char *name;
2462 PetscInt numComp;
2463
2464 ierr = PetscSectionGetFieldName(s, f, &name);CHKERRQ(ierr);
2465 ierr = PetscSectionSetFieldName(sNew, f, name);CHKERRQ(ierr);
2466 ierr = PetscSectionGetFieldComponents(s, f, &numComp);CHKERRQ(ierr);
2467 ierr = PetscSectionSetFieldComponents(sNew, f, numComp);CHKERRQ(ierr);
2468 for (c = 0; c < s->numFieldComponents[f]; ++c) {
2469 ierr = PetscSectionGetComponentName(s, f, c, &name);CHKERRQ(ierr);
2470 ierr = PetscSectionSetComponentName(sNew, f, c, name);CHKERRQ(ierr);
2471 }
2472 }
2473 ierr = ISGetLocalSize(permutation, &numPoints);CHKERRQ(ierr);
2474 ierr = ISGetIndices(permutation, &perm);CHKERRQ(ierr);
2475 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
2476 ierr = PetscSectionSetChart(sNew, pStart, pEnd);CHKERRQ(ierr);
2477 if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2478 for (p = pStart; p < pEnd; ++p) {
2479 PetscInt dof, cdof;
2480
2481 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
2482 ierr = PetscSectionSetDof(sNew, perm[p], dof);CHKERRQ(ierr);
2483 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
2484 if (cdof) {ierr = PetscSectionSetConstraintDof(sNew, perm[p], cdof);CHKERRQ(ierr);}
2485 for (f = 0; f < numFields; ++f) {
2486 ierr = PetscSectionGetFieldDof(s, p, f, &dof);CHKERRQ(ierr);
2487 ierr = PetscSectionSetFieldDof(sNew, perm[p], f, dof);CHKERRQ(ierr);
2488 ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr);
2489 if (cdof) {ierr = PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);CHKERRQ(ierr);}
2490 }
2491 }
2492 ierr = PetscSectionSetUp(sNew);CHKERRQ(ierr);
2493 for (p = pStart; p < pEnd; ++p) {
2494 const PetscInt *cind;
2495 PetscInt cdof;
2496
2497 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
2498 if (cdof) {
2499 ierr = PetscSectionGetConstraintIndices(s, p, &cind);CHKERRQ(ierr);
2500 ierr = PetscSectionSetConstraintIndices(sNew, perm[p], cind);CHKERRQ(ierr);
2501 }
2502 for (f = 0; f < numFields; ++f) {
2503 ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr);
2504 if (cdof) {
2505 ierr = PetscSectionGetFieldConstraintIndices(s, p, f, &cind);CHKERRQ(ierr);
2506 ierr = PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);CHKERRQ(ierr);
2507 }
2508 }
2509 }
2510 ierr = ISRestoreIndices(permutation, &perm);CHKERRQ(ierr);
2511 *sectionNew = sNew;
2512 PetscFunctionReturn(0);
2513 }
2514
2515 /* TODO: the next three functions should be moved to sf/utils */
2516 #include <petsc/private/sfimpl.h>
2517
2518 /*@C
2519 PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2520
2521 Collective on sf
2522
2523 Input Parameters:
2524 + sf - The SF
2525 - rootSection - Section defined on root space
2526
2527 Output Parameters:
2528 + remoteOffsets - root offsets in leaf storage, or NULL
2529 - leafSection - Section defined on the leaf space
2530
2531 Level: advanced
2532
2533 .seealso: PetscSFCreate()
2534 @*/
PetscSFDistributeSection(PetscSF sf,PetscSection rootSection,PetscInt ** remoteOffsets,PetscSection leafSection)2535 PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2536 {
2537 PetscSF embedSF;
2538 const PetscInt *indices;
2539 IS selected;
2540 PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
2541 PetscBool *sub, hasc;
2542 PetscErrorCode ierr;
2543
2544 PetscFunctionBegin;
2545 ierr = PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr);
2546 ierr = PetscSectionGetNumFields(rootSection, &numFields);CHKERRQ(ierr);
2547 if (numFields) {ierr = PetscSectionSetNumFields(leafSection, numFields);CHKERRQ(ierr);}
2548 ierr = PetscMalloc1(numFields+2, &sub);CHKERRQ(ierr);
2549 sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2550 for (f = 0; f < numFields; ++f) {
2551 PetscSectionSym sym;
2552 const char *name = NULL;
2553 PetscInt numComp = 0;
2554
2555 sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2556 ierr = PetscSectionGetFieldComponents(rootSection, f, &numComp);CHKERRQ(ierr);
2557 ierr = PetscSectionGetFieldName(rootSection, f, &name);CHKERRQ(ierr);
2558 ierr = PetscSectionGetFieldSym(rootSection, f, &sym);CHKERRQ(ierr);
2559 ierr = PetscSectionSetFieldComponents(leafSection, f, numComp);CHKERRQ(ierr);
2560 ierr = PetscSectionSetFieldName(leafSection, f, name);CHKERRQ(ierr);
2561 ierr = PetscSectionSetFieldSym(leafSection, f, sym);CHKERRQ(ierr);
2562 for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
2563 ierr = PetscSectionGetComponentName(rootSection, f, c, &name);CHKERRQ(ierr);
2564 ierr = PetscSectionSetComponentName(leafSection, f, c, name);CHKERRQ(ierr);
2565 }
2566 }
2567 ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr);
2568 ierr = PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);CHKERRQ(ierr);
2569 rpEnd = PetscMin(rpEnd,nroots);
2570 rpEnd = PetscMax(rpStart,rpEnd);
2571 /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2572 sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2573 ierr = MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));CHKERRQ(ierr);
2574 if (sub[0]) {
2575 ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr);
2576 ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr);
2577 ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr);
2578 ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr);
2579 ierr = ISDestroy(&selected);CHKERRQ(ierr);
2580 } else {
2581 ierr = PetscObjectReference((PetscObject)sf);CHKERRQ(ierr);
2582 embedSF = sf;
2583 }
2584 ierr = PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);CHKERRQ(ierr);
2585 lpEnd++;
2586
2587 ierr = PetscSectionSetChart(leafSection, lpStart, lpEnd);CHKERRQ(ierr);
2588
2589 /* Constrained dof section */
2590 hasc = sub[1];
2591 for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
2592
2593 /* Could fuse these at the cost of copies and extra allocation */
2594 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr);
2595 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr);
2596 if (sub[1]) {
2597 ierr = PetscSectionCheckConstraints_Static(rootSection);CHKERRQ(ierr);
2598 ierr = PetscSectionCheckConstraints_Static(leafSection);CHKERRQ(ierr);
2599 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
2600 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
2601 }
2602 for (f = 0; f < numFields; ++f) {
2603 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr);
2604 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr);
2605 if (sub[2+f]) {
2606 ierr = PetscSectionCheckConstraints_Static(rootSection->field[f]);CHKERRQ(ierr);
2607 ierr = PetscSectionCheckConstraints_Static(leafSection->field[f]);CHKERRQ(ierr);
2608 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
2609 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
2610 }
2611 }
2612 if (remoteOffsets) {
2613 ierr = PetscMalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr);
2614 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
2615 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
2616 }
2617 ierr = PetscSectionSetUp(leafSection);CHKERRQ(ierr);
2618 if (hasc) { /* need to communicate bcIndices */
2619 PetscSF bcSF;
2620 PetscInt *rOffBc;
2621
2622 ierr = PetscMalloc1(lpEnd - lpStart, &rOffBc);CHKERRQ(ierr);
2623 if (sub[1]) {
2624 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
2625 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
2626 ierr = PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);CHKERRQ(ierr);
2627 ierr = PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);CHKERRQ(ierr);
2628 ierr = PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);CHKERRQ(ierr);
2629 ierr = PetscSFDestroy(&bcSF);CHKERRQ(ierr);
2630 }
2631 for (f = 0; f < numFields; ++f) {
2632 if (sub[2+f]) {
2633 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
2634 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
2635 ierr = PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);CHKERRQ(ierr);
2636 ierr = PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);CHKERRQ(ierr);
2637 ierr = PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);CHKERRQ(ierr);
2638 ierr = PetscSFDestroy(&bcSF);CHKERRQ(ierr);
2639 }
2640 }
2641 ierr = PetscFree(rOffBc);CHKERRQ(ierr);
2642 }
2643 ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr);
2644 ierr = PetscFree(sub);CHKERRQ(ierr);
2645 ierr = PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr);
2646 PetscFunctionReturn(0);
2647 }
2648
2649 /*@C
2650 PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
2651
2652 Collective on sf
2653
2654 Input Parameters:
2655 + sf - The SF
2656 . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2657 - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
2658
2659 Output Parameter:
2660 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2661
2662 Level: developer
2663
2664 .seealso: PetscSFCreate()
2665 @*/
PetscSFCreateRemoteOffsets(PetscSF sf,PetscSection rootSection,PetscSection leafSection,PetscInt ** remoteOffsets)2666 PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2667 {
2668 PetscSF embedSF;
2669 const PetscInt *indices;
2670 IS selected;
2671 PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2672 PetscErrorCode ierr;
2673
2674 PetscFunctionBegin;
2675 *remoteOffsets = NULL;
2676 ierr = PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);CHKERRQ(ierr);
2677 if (numRoots < 0) PetscFunctionReturn(0);
2678 ierr = PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr);
2679 ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr);
2680 ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr);
2681 ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr);
2682 ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr);
2683 ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr);
2684 ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr);
2685 ierr = ISDestroy(&selected);CHKERRQ(ierr);
2686 ierr = PetscCalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr);
2687 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
2688 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
2689 ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr);
2690 ierr = PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr);
2691 PetscFunctionReturn(0);
2692 }
2693
2694 /*@C
2695 PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2696
2697 Collective on sf
2698
2699 Input Parameters:
2700 + sf - The SF
2701 . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2702 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2703 - leafSection - Data layout of local points for incoming data (this is the distributed section)
2704
2705 Output Parameters:
2706 - sectionSF - The new SF
2707
2708 Note: Either rootSection or remoteOffsets can be specified
2709
2710 Level: advanced
2711
2712 .seealso: PetscSFCreate()
2713 @*/
PetscSFCreateSectionSF(PetscSF sf,PetscSection rootSection,PetscInt remoteOffsets[],PetscSection leafSection,PetscSF * sectionSF)2714 PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2715 {
2716 MPI_Comm comm;
2717 const PetscInt *localPoints;
2718 const PetscSFNode *remotePoints;
2719 PetscInt lpStart, lpEnd;
2720 PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2721 PetscInt *localIndices;
2722 PetscSFNode *remoteIndices;
2723 PetscInt i, ind;
2724 PetscErrorCode ierr;
2725
2726 PetscFunctionBegin;
2727 PetscValidHeaderSpecific(sf,PETSCSF_CLASSID,1);
2728 PetscValidPointer(rootSection,2);
2729 /* Cannot check PetscValidIntPointer(remoteOffsets,3) because it can be NULL if sf does not reference any points in leafSection */
2730 PetscValidPointer(leafSection,4);
2731 PetscValidPointer(sectionSF,5);
2732 ierr = PetscObjectGetComm((PetscObject)sf,&comm);CHKERRQ(ierr);
2733 ierr = PetscSFCreate(comm, sectionSF);CHKERRQ(ierr);
2734 ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr);
2735 ierr = PetscSectionGetStorageSize(rootSection, &numSectionRoots);CHKERRQ(ierr);
2736 ierr = PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);CHKERRQ(ierr);
2737 if (numRoots < 0) PetscFunctionReturn(0);
2738 ierr = PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr);
2739 for (i = 0; i < numPoints; ++i) {
2740 PetscInt localPoint = localPoints ? localPoints[i] : i;
2741 PetscInt dof;
2742
2743 if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2744 ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr);
2745 numIndices += dof;
2746 }
2747 }
2748 ierr = PetscMalloc1(numIndices, &localIndices);CHKERRQ(ierr);
2749 ierr = PetscMalloc1(numIndices, &remoteIndices);CHKERRQ(ierr);
2750 /* Create new index graph */
2751 for (i = 0, ind = 0; i < numPoints; ++i) {
2752 PetscInt localPoint = localPoints ? localPoints[i] : i;
2753 PetscInt rank = remotePoints[i].rank;
2754
2755 if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2756 PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2757 PetscInt loff, dof, d;
2758
2759 ierr = PetscSectionGetOffset(leafSection, localPoint, &loff);CHKERRQ(ierr);
2760 ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr);
2761 for (d = 0; d < dof; ++d, ++ind) {
2762 localIndices[ind] = loff+d;
2763 remoteIndices[ind].rank = rank;
2764 remoteIndices[ind].index = remoteOffset+d;
2765 }
2766 }
2767 }
2768 if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2769 ierr = PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);CHKERRQ(ierr);
2770 ierr = PetscSFSetUp(*sectionSF);CHKERRQ(ierr);
2771 ierr = PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr);
2772 PetscFunctionReturn(0);
2773 }
2774
2775 /*@
2776 PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2777
2778 Collective on section
2779
2780 Input Parameters:
2781 + section - The PetscSection
2782 . obj - A PetscObject which serves as the key for this index
2783 . clSection - Section giving the size of the closure of each point
2784 - clPoints - IS giving the points in each closure
2785
2786 Note: We compress out closure points with no dofs in this section
2787
2788 Level: advanced
2789
2790 .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2791 @*/
PetscSectionSetClosureIndex(PetscSection section,PetscObject obj,PetscSection clSection,IS clPoints)2792 PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2793 {
2794 PetscErrorCode ierr;
2795
2796 PetscFunctionBegin;
2797 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
2798 PetscValidHeaderSpecific(clSection,PETSC_SECTION_CLASSID,3);
2799 PetscValidHeaderSpecific(clPoints,IS_CLASSID,4);
2800 if (section->clObj != obj) {ierr = PetscSectionResetClosurePermutation(section);CHKERRQ(ierr);}
2801 section->clObj = obj;
2802 ierr = PetscObjectReference((PetscObject)clSection);CHKERRQ(ierr);
2803 ierr = PetscObjectReference((PetscObject)clPoints);CHKERRQ(ierr);
2804 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr);
2805 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr);
2806 section->clSection = clSection;
2807 section->clPoints = clPoints;
2808 PetscFunctionReturn(0);
2809 }
2810
2811 /*@
2812 PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2813
2814 Collective on section
2815
2816 Input Parameters:
2817 + section - The PetscSection
2818 - obj - A PetscObject which serves as the key for this index
2819
2820 Output Parameters:
2821 + clSection - Section giving the size of the closure of each point
2822 - clPoints - IS giving the points in each closure
2823
2824 Note: We compress out closure points with no dofs in this section
2825
2826 Level: advanced
2827
2828 .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2829 @*/
PetscSectionGetClosureIndex(PetscSection section,PetscObject obj,PetscSection * clSection,IS * clPoints)2830 PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2831 {
2832 PetscFunctionBegin;
2833 if (section->clObj == obj) {
2834 if (clSection) *clSection = section->clSection;
2835 if (clPoints) *clPoints = section->clPoints;
2836 } else {
2837 if (clSection) *clSection = NULL;
2838 if (clPoints) *clPoints = NULL;
2839 }
2840 PetscFunctionReturn(0);
2841 }
2842
PetscSectionSetClosurePermutation_Internal(PetscSection section,PetscObject obj,PetscInt depth,PetscInt clSize,PetscCopyMode mode,PetscInt * clPerm)2843 PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2844 {
2845 PetscInt i;
2846 khiter_t iter;
2847 int new_entry;
2848 PetscSectionClosurePermKey key = {depth, clSize};
2849 PetscSectionClosurePermVal *val;
2850 PetscErrorCode ierr;
2851
2852 PetscFunctionBegin;
2853 if (section->clObj != obj) {
2854 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr);
2855 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr);
2856 }
2857 section->clObj = obj;
2858 if (!section->clHash) {ierr = PetscClPermCreate(§ion->clHash);CHKERRQ(ierr);}
2859 iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2860 val = &kh_val(section->clHash, iter);
2861 if (!new_entry) {
2862 ierr = PetscFree(val->perm);CHKERRQ(ierr);
2863 ierr = PetscFree(val->invPerm);CHKERRQ(ierr);
2864 }
2865 if (mode == PETSC_COPY_VALUES) {
2866 ierr = PetscMalloc1(clSize, &val->perm);CHKERRQ(ierr);
2867 ierr = PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));CHKERRQ(ierr);
2868 ierr = PetscArraycpy(val->perm, clPerm, clSize);CHKERRQ(ierr);
2869 } else if (mode == PETSC_OWN_POINTER) {
2870 val->perm = clPerm;
2871 } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2872 ierr = PetscMalloc1(clSize, &val->invPerm);CHKERRQ(ierr);
2873 for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2874 PetscFunctionReturn(0);
2875 }
2876
2877 /*@
2878 PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2879
2880 Not Collective
2881
2882 Input Parameters:
2883 + section - The PetscSection
2884 . obj - A PetscObject which serves as the key for this index (usually a DM)
2885 . depth - Depth of points on which to apply the given permutation
2886 - perm - Permutation of the cell dof closure
2887
2888 Note:
2889 The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2890 mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2891 each topology and degree.
2892
2893 This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2894 exotic/enriched spaces on mixed topology meshes.
2895
2896 Level: intermediate
2897
2898 .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2899 @*/
PetscSectionSetClosurePermutation(PetscSection section,PetscObject obj,PetscInt depth,IS perm)2900 PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2901 {
2902 const PetscInt *clPerm = NULL;
2903 PetscInt clSize = 0;
2904 PetscErrorCode ierr;
2905
2906 PetscFunctionBegin;
2907 if (perm) {
2908 ierr = ISGetLocalSize(perm, &clSize);CHKERRQ(ierr);
2909 ierr = ISGetIndices(perm, &clPerm);CHKERRQ(ierr);
2910 }
2911 ierr = PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);CHKERRQ(ierr);
2912 if (perm) {ierr = ISRestoreIndices(perm, &clPerm);CHKERRQ(ierr);}
2913 PetscFunctionReturn(0);
2914 }
2915
PetscSectionGetClosurePermutation_Internal(PetscSection section,PetscObject obj,PetscInt depth,PetscInt size,const PetscInt * perm[])2916 PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2917 {
2918 PetscErrorCode ierr;
2919
2920 PetscFunctionBegin;
2921 if (section->clObj == obj) {
2922 PetscSectionClosurePermKey k = {depth, size};
2923 PetscSectionClosurePermVal v;
2924 ierr = PetscClPermGet(section->clHash, k, &v);CHKERRQ(ierr);
2925 if (perm) *perm = v.perm;
2926 } else {
2927 if (perm) *perm = NULL;
2928 }
2929 PetscFunctionReturn(0);
2930 }
2931
2932 /*@
2933 PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2934
2935 Not collective
2936
2937 Input Parameters:
2938 + section - The PetscSection
2939 . obj - A PetscObject which serves as the key for this index (usually a DM)
2940 . depth - Depth stratum on which to obtain closure permutation
2941 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2942
2943 Output Parameter:
2944 . perm - The dof closure permutation
2945
2946 Note:
2947 The user must destroy the IS that is returned.
2948
2949 Level: intermediate
2950
2951 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2952 @*/
PetscSectionGetClosurePermutation(PetscSection section,PetscObject obj,PetscInt depth,PetscInt clSize,IS * perm)2953 PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2954 {
2955 const PetscInt *clPerm;
2956 PetscErrorCode ierr;
2957
2958 PetscFunctionBegin;
2959 ierr = PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);CHKERRQ(ierr);
2960 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr);
2961 PetscFunctionReturn(0);
2962 }
2963
PetscSectionGetClosureInversePermutation_Internal(PetscSection section,PetscObject obj,PetscInt depth,PetscInt size,const PetscInt * perm[])2964 PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2965 {
2966 PetscErrorCode ierr;
2967
2968 PetscFunctionBegin;
2969 if (section->clObj == obj && section->clHash) {
2970 PetscSectionClosurePermKey k = {depth, size};
2971 PetscSectionClosurePermVal v;
2972 ierr = PetscClPermGet(section->clHash, k, &v);CHKERRQ(ierr);
2973 if (perm) *perm = v.invPerm;
2974 } else {
2975 if (perm) *perm = NULL;
2976 }
2977 PetscFunctionReturn(0);
2978 }
2979
2980 /*@
2981 PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2982
2983 Not collective
2984
2985 Input Parameters:
2986 + section - The PetscSection
2987 . obj - A PetscObject which serves as the key for this index (usually a DM)
2988 . depth - Depth stratum on which to obtain closure permutation
2989 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2990
2991 Output Parameters:
2992 . perm - The dof closure permutation
2993
2994 Note:
2995 The user must destroy the IS that is returned.
2996
2997 Level: intermediate
2998
2999 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
3000 @*/
PetscSectionGetClosureInversePermutation(PetscSection section,PetscObject obj,PetscInt depth,PetscInt clSize,IS * perm)3001 PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3002 {
3003 const PetscInt *clPerm;
3004 PetscErrorCode ierr;
3005
3006 PetscFunctionBegin;
3007 ierr = PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);CHKERRQ(ierr);
3008 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr);
3009 PetscFunctionReturn(0);
3010 }
3011
3012 /*@
3013 PetscSectionGetField - Get the subsection associated with a single field
3014
3015 Input Parameters:
3016 + s - The PetscSection
3017 - field - The field number
3018
3019 Output Parameter:
3020 . subs - The subsection for the given field
3021
3022 Level: intermediate
3023
3024 .seealso: PetscSectionSetNumFields()
3025 @*/
PetscSectionGetField(PetscSection s,PetscInt field,PetscSection * subs)3026 PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3027 {
3028 PetscFunctionBegin;
3029 PetscValidHeaderSpecific(s,PETSC_SECTION_CLASSID,1);
3030 PetscValidPointer(subs,3);
3031 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);
3032 *subs = s->field[field];
3033 PetscFunctionReturn(0);
3034 }
3035
3036 PetscClassId PETSC_SECTION_SYM_CLASSID;
3037 PetscFunctionList PetscSectionSymList = NULL;
3038
3039 /*@
3040 PetscSectionSymCreate - Creates an empty PetscSectionSym object.
3041
3042 Collective
3043
3044 Input Parameter:
3045 . comm - the MPI communicator
3046
3047 Output Parameter:
3048 . sym - pointer to the new set of symmetries
3049
3050 Level: developer
3051
3052 .seealso: PetscSectionSym, PetscSectionSymDestroy()
3053 @*/
PetscSectionSymCreate(MPI_Comm comm,PetscSectionSym * sym)3054 PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3055 {
3056 PetscErrorCode ierr;
3057
3058 PetscFunctionBegin;
3059 PetscValidPointer(sym,2);
3060 ierr = ISInitializePackage();CHKERRQ(ierr);
3061 ierr = PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);CHKERRQ(ierr);
3062 PetscFunctionReturn(0);
3063 }
3064
3065 /*@C
3066 PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
3067
3068 Collective on PetscSectionSym
3069
3070 Input Parameters:
3071 + sym - The section symmetry object
3072 - method - The name of the section symmetry type
3073
3074 Level: developer
3075
3076 .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
3077 @*/
PetscSectionSymSetType(PetscSectionSym sym,PetscSectionSymType method)3078 PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3079 {
3080 PetscErrorCode (*r)(PetscSectionSym);
3081 PetscBool match;
3082 PetscErrorCode ierr;
3083
3084 PetscFunctionBegin;
3085 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1);
3086 ierr = PetscObjectTypeCompare((PetscObject) sym, method, &match);CHKERRQ(ierr);
3087 if (match) PetscFunctionReturn(0);
3088
3089 ierr = PetscFunctionListFind(PetscSectionSymList,method,&r);CHKERRQ(ierr);
3090 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3091 if (sym->ops->destroy) {
3092 ierr = (*sym->ops->destroy)(sym);CHKERRQ(ierr);
3093 sym->ops->destroy = NULL;
3094 }
3095 ierr = (*r)(sym);CHKERRQ(ierr);
3096 ierr = PetscObjectChangeTypeName((PetscObject)sym,method);CHKERRQ(ierr);
3097 PetscFunctionReturn(0);
3098 }
3099
3100
3101 /*@C
3102 PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
3103
3104 Not Collective
3105
3106 Input Parameter:
3107 . sym - The section symmetry
3108
3109 Output Parameter:
3110 . type - The index set type name
3111
3112 Level: developer
3113
3114 .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
3115 @*/
PetscSectionSymGetType(PetscSectionSym sym,PetscSectionSymType * type)3116 PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3117 {
3118 PetscFunctionBegin;
3119 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1);
3120 PetscValidCharPointer(type,2);
3121 *type = ((PetscObject)sym)->type_name;
3122 PetscFunctionReturn(0);
3123 }
3124
3125 /*@C
3126 PetscSectionSymRegister - Adds a new section symmetry implementation
3127
3128 Not Collective
3129
3130 Input Parameters:
3131 + name - The name of a new user-defined creation routine
3132 - create_func - The creation routine itself
3133
3134 Notes:
3135 PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
3136
3137 Level: developer
3138
3139 .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3140 @*/
PetscSectionSymRegister(const char sname[],PetscErrorCode (* function)(PetscSectionSym))3141 PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3142 {
3143 PetscErrorCode ierr;
3144
3145 PetscFunctionBegin;
3146 ierr = ISInitializePackage();CHKERRQ(ierr);
3147 ierr = PetscFunctionListAdd(&PetscSectionSymList,sname,function);CHKERRQ(ierr);
3148 PetscFunctionReturn(0);
3149 }
3150
3151 /*@
3152 PetscSectionSymDestroy - Destroys a section symmetry.
3153
3154 Collective on PetscSectionSym
3155
3156 Input Parameters:
3157 . sym - the section symmetry
3158
3159 Level: developer
3160
3161 .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3162 @*/
PetscSectionSymDestroy(PetscSectionSym * sym)3163 PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3164 {
3165 SymWorkLink link,next;
3166 PetscErrorCode ierr;
3167
3168 PetscFunctionBegin;
3169 if (!*sym) PetscFunctionReturn(0);
3170 PetscValidHeaderSpecific((*sym),PETSC_SECTION_SYM_CLASSID,1);
3171 if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; PetscFunctionReturn(0);}
3172 if ((*sym)->ops->destroy) {
3173 ierr = (*(*sym)->ops->destroy)(*sym);CHKERRQ(ierr);
3174 }
3175 if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3176 for (link=(*sym)->workin; link; link=next) {
3177 next = link->next;
3178 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);CHKERRQ(ierr);
3179 ierr = PetscFree(link);CHKERRQ(ierr);
3180 }
3181 (*sym)->workin = NULL;
3182 ierr = PetscHeaderDestroy(sym);CHKERRQ(ierr);
3183 PetscFunctionReturn(0);
3184 }
3185
3186 /*@C
3187 PetscSectionSymView - Displays a section symmetry
3188
3189 Collective on PetscSectionSym
3190
3191 Input Parameters:
3192 + sym - the index set
3193 - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
3194
3195 Level: developer
3196
3197 .seealso: PetscViewerASCIIOpen()
3198 @*/
PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)3199 PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3200 {
3201 PetscErrorCode ierr;
3202
3203 PetscFunctionBegin;
3204 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1);
3205 if (!viewer) {
3206 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);CHKERRQ(ierr);
3207 }
3208 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
3209 PetscCheckSameComm(sym,1,viewer,2);
3210 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);CHKERRQ(ierr);
3211 if (sym->ops->view) {
3212 ierr = (*sym->ops->view)(sym,viewer);CHKERRQ(ierr);
3213 }
3214 PetscFunctionReturn(0);
3215 }
3216
3217 /*@
3218 PetscSectionSetSym - Set the symmetries for the data referred to by the section
3219
3220 Collective on PetscSection
3221
3222 Input Parameters:
3223 + section - the section describing data layout
3224 - sym - the symmetry describing the affect of orientation on the access of the data
3225
3226 Level: developer
3227
3228 .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3229 @*/
PetscSectionSetSym(PetscSection section,PetscSectionSym sym)3230 PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3231 {
3232 PetscErrorCode ierr;
3233
3234 PetscFunctionBegin;
3235 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3236 ierr = PetscSectionSymDestroy(&(section->sym));CHKERRQ(ierr);
3237 if (sym) {
3238 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,2);
3239 PetscCheckSameComm(section,1,sym,2);
3240 ierr = PetscObjectReference((PetscObject) sym);CHKERRQ(ierr);
3241 }
3242 section->sym = sym;
3243 PetscFunctionReturn(0);
3244 }
3245
3246 /*@
3247 PetscSectionGetSym - Get the symmetries for the data referred to by the section
3248
3249 Not collective
3250
3251 Input Parameters:
3252 . section - the section describing data layout
3253
3254 Output Parameters:
3255 . sym - the symmetry describing the affect of orientation on the access of the data
3256
3257 Level: developer
3258
3259 .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3260 @*/
PetscSectionGetSym(PetscSection section,PetscSectionSym * sym)3261 PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3262 {
3263 PetscFunctionBegin;
3264 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3265 *sym = section->sym;
3266 PetscFunctionReturn(0);
3267 }
3268
3269 /*@
3270 PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3271
3272 Collective on PetscSection
3273
3274 Input Parameters:
3275 + section - the section describing data layout
3276 . field - the field number
3277 - sym - the symmetry describing the affect of orientation on the access of the data
3278
3279 Level: developer
3280
3281 .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3282 @*/
PetscSectionSetFieldSym(PetscSection section,PetscInt field,PetscSectionSym sym)3283 PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3284 {
3285 PetscErrorCode ierr;
3286
3287 PetscFunctionBegin;
3288 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3289 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);
3290 ierr = PetscSectionSetSym(section->field[field],sym);CHKERRQ(ierr);
3291 PetscFunctionReturn(0);
3292 }
3293
3294 /*@
3295 PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3296
3297 Collective on PetscSection
3298
3299 Input Parameters:
3300 + section - the section describing data layout
3301 - field - the field number
3302
3303 Output Parameters:
3304 . sym - the symmetry describing the affect of orientation on the access of the data
3305
3306 Level: developer
3307
3308 .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3309 @*/
PetscSectionGetFieldSym(PetscSection section,PetscInt field,PetscSectionSym * sym)3310 PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3311 {
3312 PetscFunctionBegin;
3313 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3314 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);
3315 *sym = section->field[field]->sym;
3316 PetscFunctionReturn(0);
3317 }
3318
3319 /*@C
3320 PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3321
3322 Not collective
3323
3324 Input Parameters:
3325 + section - the section
3326 . numPoints - the number of points
3327 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3328 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3329 context, see DMPlexGetConeOrientation()).
3330
3331 Output Parameter:
3332 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3333 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3334 identity).
3335
3336 Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3337 .vb
3338 const PetscInt **perms;
3339 const PetscScalar **rots;
3340 PetscInt lOffset;
3341
3342 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3343 for (i = 0, lOffset = 0; i < numPoints; i++) {
3344 PetscInt point = points[2*i], dof, sOffset;
3345 const PetscInt *perm = perms ? perms[i] : NULL;
3346 const PetscScalar *rot = rots ? rots[i] : NULL;
3347
3348 PetscSectionGetDof(section,point,&dof);
3349 PetscSectionGetOffset(section,point,&sOffset);
3350
3351 if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3352 else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3353 if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3354 lOffset += dof;
3355 }
3356 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3357 .ve
3358
3359 Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3360 .vb
3361 const PetscInt **perms;
3362 const PetscScalar **rots;
3363 PetscInt lOffset;
3364
3365 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3366 for (i = 0, lOffset = 0; i < numPoints; i++) {
3367 PetscInt point = points[2*i], dof, sOffset;
3368 const PetscInt *perm = perms ? perms[i] : NULL;
3369 const PetscScalar *rot = rots ? rots[i] : NULL;
3370
3371 PetscSectionGetDof(section,point,&dof);
3372 PetscSectionGetOffset(section,point,&sOff);
3373
3374 if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3375 else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3376 offset += dof;
3377 }
3378 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3379 .ve
3380
3381 Level: developer
3382
3383 .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3384 @*/
PetscSectionGetPointSyms(PetscSection section,PetscInt numPoints,const PetscInt * points,const PetscInt *** perms,const PetscScalar *** rots)3385 PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3386 {
3387 PetscSectionSym sym;
3388 PetscErrorCode ierr;
3389
3390 PetscFunctionBegin;
3391 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3392 if (numPoints) PetscValidIntPointer(points,3);
3393 if (perms) *perms = NULL;
3394 if (rots) *rots = NULL;
3395 sym = section->sym;
3396 if (sym && (perms || rots)) {
3397 SymWorkLink link;
3398
3399 if (sym->workin) {
3400 link = sym->workin;
3401 sym->workin = sym->workin->next;
3402 } else {
3403 ierr = PetscNewLog(sym,&link);CHKERRQ(ierr);
3404 }
3405 if (numPoints > link->numPoints) {
3406 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);CHKERRQ(ierr);
3407 ierr = PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);CHKERRQ(ierr);
3408 link->numPoints = numPoints;
3409 }
3410 link->next = sym->workout;
3411 sym->workout = link;
3412 ierr = PetscArrayzero((PetscInt**)link->perms,numPoints);CHKERRQ(ierr);
3413 ierr = PetscArrayzero((PetscInt**)link->rots,numPoints);CHKERRQ(ierr);
3414 ierr = (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);CHKERRQ(ierr);
3415 if (perms) *perms = link->perms;
3416 if (rots) *rots = link->rots;
3417 }
3418 PetscFunctionReturn(0);
3419 }
3420
3421 /*@C
3422 PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3423
3424 Not collective
3425
3426 Input Parameters:
3427 + section - the section
3428 . numPoints - the number of points
3429 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3430 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3431 context, see DMPlexGetConeOrientation()).
3432
3433 Output Parameter:
3434 + perms - The permutations for the given orientations: set to NULL at conclusion
3435 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3436
3437 Level: developer
3438
3439 .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3440 @*/
PetscSectionRestorePointSyms(PetscSection section,PetscInt numPoints,const PetscInt * points,const PetscInt *** perms,const PetscScalar *** rots)3441 PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3442 {
3443 PetscSectionSym sym;
3444
3445 PetscFunctionBegin;
3446 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3447 sym = section->sym;
3448 if (sym && (perms || rots)) {
3449 SymWorkLink *p,link;
3450
3451 for (p=&sym->workout; (link=*p); p=&link->next) {
3452 if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3453 *p = link->next;
3454 link->next = sym->workin;
3455 sym->workin = link;
3456 if (perms) *perms = NULL;
3457 if (rots) *rots = NULL;
3458 PetscFunctionReturn(0);
3459 }
3460 }
3461 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3462 }
3463 PetscFunctionReturn(0);
3464 }
3465
3466 /*@C
3467 PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3468
3469 Not collective
3470
3471 Input Parameters:
3472 + section - the section
3473 . field - the field of the section
3474 . numPoints - the number of points
3475 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3476 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3477 context, see DMPlexGetConeOrientation()).
3478
3479 Output Parameter:
3480 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3481 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3482 identity).
3483
3484 Level: developer
3485
3486 .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3487 @*/
PetscSectionGetFieldPointSyms(PetscSection section,PetscInt field,PetscInt numPoints,const PetscInt * points,const PetscInt *** perms,const PetscScalar *** rots)3488 PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3489 {
3490 PetscErrorCode ierr;
3491
3492 PetscFunctionBegin;
3493 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3494 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);
3495 ierr = PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr);
3496 PetscFunctionReturn(0);
3497 }
3498
3499 /*@C
3500 PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3501
3502 Not collective
3503
3504 Input Parameters:
3505 + section - the section
3506 . field - the field number
3507 . numPoints - the number of points
3508 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3509 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3510 context, see DMPlexGetConeOrientation()).
3511
3512 Output Parameter:
3513 + perms - The permutations for the given orientations: set to NULL at conclusion
3514 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3515
3516 Level: developer
3517
3518 .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3519 @*/
PetscSectionRestoreFieldPointSyms(PetscSection section,PetscInt field,PetscInt numPoints,const PetscInt * points,const PetscInt *** perms,const PetscScalar *** rots)3520 PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3521 {
3522 PetscErrorCode ierr;
3523
3524 PetscFunctionBegin;
3525 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1);
3526 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);
3527 ierr = PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr);
3528 PetscFunctionReturn(0);
3529 }
3530
3531 /*@
3532 PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3533
3534 Not collective
3535
3536 Input Parameter:
3537 . s - the global PetscSection
3538
3539 Output Parameters:
3540 . flg - the flag
3541
3542 Level: developer
3543
3544 .seealso: PetscSectionSetChart(), PetscSectionCreate()
3545 @*/
PetscSectionGetUseFieldOffsets(PetscSection s,PetscBool * flg)3546 PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3547 {
3548 PetscFunctionBegin;
3549 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
3550 *flg = s->useFieldOff;
3551 PetscFunctionReturn(0);
3552 }
3553
3554 /*@
3555 PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3556
3557 Not collective
3558
3559 Input Parameters:
3560 + s - the global PetscSection
3561 - flg - the flag
3562
3563 Level: developer
3564
3565 .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3566 @*/
PetscSectionSetUseFieldOffsets(PetscSection s,PetscBool flg)3567 PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3568 {
3569 PetscFunctionBegin;
3570 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
3571 s->useFieldOff = flg;
3572 PetscFunctionReturn(0);
3573 }
3574
3575 #define PetscSectionExpandPoints_Loop(TYPE) \
3576 { \
3577 PetscInt i, n, o0, o1, size; \
3578 TYPE *a0 = (TYPE*)origArray, *a1; \
3579 ierr = PetscSectionGetStorageSize(s, &size);CHKERRQ(ierr); \
3580 ierr = PetscMalloc1(size, &a1);CHKERRQ(ierr); \
3581 for (i=0; i<npoints; i++) { \
3582 ierr = PetscSectionGetOffset(origSection, points_[i], &o0);CHKERRQ(ierr); \
3583 ierr = PetscSectionGetOffset(s, i, &o1);CHKERRQ(ierr); \
3584 ierr = PetscSectionGetDof(s, i, &n);CHKERRQ(ierr); \
3585 ierr = PetscMemcpy(&a1[o1], &a0[o0], n*unitsize);CHKERRQ(ierr); \
3586 } \
3587 *newArray = (void*)a1; \
3588 }
3589
3590 /*@
3591 PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3592
3593 Not collective
3594
3595 Input Parameters:
3596 + origSection - the PetscSection describing the layout of the array
3597 . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3598 . origArray - the array; its size must be equal to the storage size of origSection
3599 - points - IS with points to extract; its indices must lie in the chart of origSection
3600
3601 Output Parameters:
3602 + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3603 - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3604
3605 Level: developer
3606
3607 .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3608 @*/
PetscSectionExtractDofsFromArray(PetscSection origSection,MPI_Datatype dataType,const void * origArray,IS points,PetscSection * newSection,void * newArray[])3609 PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3610 {
3611 PetscSection s;
3612 const PetscInt *points_;
3613 PetscInt i, n, npoints, pStart, pEnd;
3614 PetscMPIInt unitsize;
3615 PetscErrorCode ierr;
3616
3617 PetscFunctionBegin;
3618 PetscValidHeaderSpecific(origSection, PETSC_SECTION_CLASSID, 1);
3619 PetscValidPointer(origArray, 3);
3620 PetscValidHeaderSpecific(points, IS_CLASSID, 4);
3621 if (newSection) PetscValidPointer(newSection, 5);
3622 if (newArray) PetscValidPointer(newArray, 6);
3623 ierr = MPI_Type_size(dataType, &unitsize);CHKERRQ(ierr);
3624 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr);
3625 ierr = ISGetIndices(points, &points_);CHKERRQ(ierr);
3626 ierr = PetscSectionGetChart(origSection, &pStart, &pEnd);CHKERRQ(ierr);
3627 ierr = PetscSectionCreate(PETSC_COMM_SELF, &s);CHKERRQ(ierr);
3628 ierr = PetscSectionSetChart(s, 0, npoints);CHKERRQ(ierr);
3629 for (i=0; i<npoints; i++) {
3630 if (PetscUnlikely(points_[i] < pStart || points_[i] >= pEnd)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %d (index %d) in input IS out of input section's chart", points_[i], i);
3631 ierr = PetscSectionGetDof(origSection, points_[i], &n);CHKERRQ(ierr);
3632 ierr = PetscSectionSetDof(s, i, n);CHKERRQ(ierr);
3633 }
3634 ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
3635 if (newArray) {
3636 if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3637 else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3638 else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3639 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3640 }
3641 if (newSection) {
3642 *newSection = s;
3643 } else {
3644 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
3645 }
3646 PetscFunctionReturn(0);
3647 }
3648