1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/isimpl.h>
3 #include <petsc/private/petscfeimpl.h>
4 #include <petscsf.h>
5 #include <petscds.h>
6
7 /** hierarchy routines */
8
9 /*@
10 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes.
11
12 Not collective
13
14 Input Parameters:
15 + dm - The DMPlex object
16 - ref - The reference tree DMPlex object
17
18 Level: intermediate
19
20 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree()
21 @*/
DMPlexSetReferenceTree(DM dm,DM ref)22 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
23 {
24 DM_Plex *mesh = (DM_Plex *)dm->data;
25 PetscErrorCode ierr;
26
27 PetscFunctionBegin;
28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
29 if (ref) {PetscValidHeaderSpecific(ref, DM_CLASSID, 2);}
30 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr);
31 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr);
32 mesh->referenceTree = ref;
33 PetscFunctionReturn(0);
34 }
35
36 /*@
37 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.
38
39 Not collective
40
41 Input Parameters:
42 . dm - The DMPlex object
43
44 Output Parameters:
45 . ref - The reference tree DMPlex object
46
47 Level: intermediate
48
49 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree()
50 @*/
DMPlexGetReferenceTree(DM dm,DM * ref)51 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
52 {
53 DM_Plex *mesh = (DM_Plex *)dm->data;
54
55 PetscFunctionBegin;
56 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
57 PetscValidPointer(ref,2);
58 *ref = mesh->referenceTree;
59 PetscFunctionReturn(0);
60 }
61
DMPlexReferenceTreeGetChildSymmetry_Default(DM dm,PetscInt parent,PetscInt parentOrientA,PetscInt childOrientA,PetscInt childA,PetscInt parentOrientB,PetscInt * childOrientB,PetscInt * childB)62 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
63 {
64 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
65 PetscErrorCode ierr;
66
67 PetscFunctionBegin;
68 if (parentOrientA == parentOrientB) {
69 if (childOrientB) *childOrientB = childOrientA;
70 if (childB) *childB = childA;
71 PetscFunctionReturn(0);
72 }
73 for (dim = 0; dim < 3; dim++) {
74 ierr = DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);CHKERRQ(ierr);
75 if (parent >= dStart && parent <= dEnd) {
76 break;
77 }
78 }
79 if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim);
80 if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children");
81 if (childA < dStart || childA >= dEnd) {
82 /* this is a lower-dimensional child: bootstrap */
83 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
84 const PetscInt *supp, *coneA, *coneB, *oA, *oB;
85
86 ierr = DMPlexGetSupportSize(dm,childA,&size);CHKERRQ(ierr);
87 ierr = DMPlexGetSupport(dm,childA,&supp);CHKERRQ(ierr);
88
89 /* find a point sA in supp(childA) that has the same parent */
90 for (i = 0; i < size; i++) {
91 PetscInt sParent;
92
93 sA = supp[i];
94 if (sA == parent) continue;
95 ierr = DMPlexGetTreeParent(dm,sA,&sParent,NULL);CHKERRQ(ierr);
96 if (sParent == parent) {
97 break;
98 }
99 }
100 if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children");
101 /* find out which point sB is in an equivalent position to sA under
102 * parentOrientB */
103 ierr = DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);CHKERRQ(ierr);
104 ierr = DMPlexGetConeSize(dm,sA,&sConeSize);CHKERRQ(ierr);
105 ierr = DMPlexGetCone(dm,sA,&coneA);CHKERRQ(ierr);
106 ierr = DMPlexGetCone(dm,sB,&coneB);CHKERRQ(ierr);
107 ierr = DMPlexGetConeOrientation(dm,sA,&oA);CHKERRQ(ierr);
108 ierr = DMPlexGetConeOrientation(dm,sB,&oB);CHKERRQ(ierr);
109 /* step through the cone of sA in natural order */
110 for (i = 0; i < sConeSize; i++) {
111 if (coneA[i] == childA) {
112 /* if childA is at position i in coneA,
113 * then we want the point that is at sOrientB*i in coneB */
114 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize);
115 if (childB) *childB = coneB[j];
116 if (childOrientB) {
117 PetscInt oBtrue;
118
119 ierr = DMPlexGetConeSize(dm,childA,&coneSize);CHKERRQ(ierr);
120 /* compose sOrientB and oB[j] */
121 if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge");
122 /* we may have to flip an edge */
123 oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0;
124 ABswap = DihedralSwap(coneSize,oA[i],oBtrue);
125 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
126 }
127 break;
128 }
129 }
130 if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch");
131 PetscFunctionReturn(0);
132 }
133 /* get the cone size and symmetry swap */
134 ierr = DMPlexGetConeSize(dm,parent,&coneSize);CHKERRQ(ierr);
135 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
136 if (dim == 2) {
137 /* orientations refer to cones: we want them to refer to vertices:
138 * if it's a rotation, they are the same, but if the order is reversed, a
139 * permutation that puts side i first does *not* put vertex i first */
140 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
141 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
142 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
143 } else {
144 ABswapVert = ABswap;
145 }
146 if (childB) {
147 /* assume that each child corresponds to a vertex, in the same order */
148 PetscInt p, posA = -1, numChildren, i;
149 const PetscInt *children;
150
151 /* count which position the child is in */
152 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr);
153 for (i = 0; i < numChildren; i++) {
154 p = children[i];
155 if (p == childA) {
156 posA = i;
157 break;
158 }
159 }
160 if (posA >= coneSize) {
161 /* this is the triangle in the middle of a uniformly refined triangle: it
162 * is invariant */
163 if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else");
164 *childB = childA;
165 }
166 else {
167 /* figure out position B by applying ABswapVert */
168 PetscInt posB;
169
170 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize);
171 if (childB) *childB = children[posB];
172 }
173 }
174 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
175 PetscFunctionReturn(0);
176 }
177
178 /*@
179 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
180
181 Input Parameters:
182 + dm - the reference tree DMPlex object
183 . parent - the parent point
184 . parentOrientA - the reference orientation for describing the parent
185 . childOrientA - the reference orientation for describing the child
186 . childA - the reference childID for describing the child
187 - parentOrientB - the new orientation for describing the parent
188
189 Output Parameters:
190 + childOrientB - if not NULL, set to the new oreintation for describing the child
191 - childB - if not NULL, the new childID for describing the child
192
193 Level: developer
194
195 .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree()
196 @*/
DMPlexReferenceTreeGetChildSymmetry(DM dm,PetscInt parent,PetscInt parentOrientA,PetscInt childOrientA,PetscInt childA,PetscInt parentOrientB,PetscInt * childOrientB,PetscInt * childB)197 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
198 {
199 DM_Plex *mesh = (DM_Plex *)dm->data;
200 PetscErrorCode ierr;
201
202 PetscFunctionBegin;
203 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
204 if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
205 ierr = mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);CHKERRQ(ierr);
206 PetscFunctionReturn(0);
207 }
208
209 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool);
210
DMPlexCreateReferenceTree_SetTree(DM dm,PetscSection parentSection,PetscInt parents[],PetscInt childIDs[])211 PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
212 {
213 PetscErrorCode ierr;
214
215 PetscFunctionBegin;
216 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
217 PetscFunctionReturn(0);
218 }
219
DMPlexCreateReferenceTree_Union(DM K,DM Kref,const char * labelName,DM * ref)220 PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
221 {
222 MPI_Comm comm;
223 PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
224 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
225 DMLabel identity, identityRef;
226 PetscSection unionSection, unionConeSection, parentSection;
227 PetscScalar *unionCoords;
228 IS perm;
229 PetscErrorCode ierr;
230
231 PetscFunctionBegin;
232 comm = PetscObjectComm((PetscObject)K);
233 ierr = DMGetDimension(K, &dim);CHKERRQ(ierr);
234 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr);
235 ierr = DMGetLabel(K, labelName, &identity);CHKERRQ(ierr);
236 ierr = DMGetLabel(Kref, labelName, &identityRef);CHKERRQ(ierr);
237 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr);
238 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr);
239 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr);
240 /* count points that will go in the union */
241 for (p = pStart; p < pEnd; p++) {
242 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr);
243 }
244 for (p = pRefStart; p < pRefEnd; p++) {
245 PetscInt q, qSize;
246 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr);
247 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr);
248 if (qSize > 1) {
249 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr);
250 }
251 }
252 ierr = PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);CHKERRQ(ierr);
253 offset = 0;
254 /* stratify points in the union by topological dimension */
255 for (d = 0; d <= dim; d++) {
256 PetscInt cStart, cEnd, c;
257
258 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr);
259 for (c = cStart; c < cEnd; c++) {
260 permvals[offset++] = c;
261 }
262
263 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr);
264 for (c = cStart; c < cEnd; c++) {
265 permvals[offset++] = c + (pEnd - pStart);
266 }
267 }
268 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr);
269 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr);
270 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr);
271 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr);
272 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr);
273 /* count dimension points */
274 for (d = 0; d <= dim; d++) {
275 PetscInt cStart, cOff, cOff2;
276 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr);
277 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr);
278 if (d < dim) {
279 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr);
280 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr);
281 }
282 else {
283 cOff2 = numUnionPoints;
284 }
285 numDimPoints[dim - d] = cOff2 - cOff;
286 }
287 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr);
288 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr);
289 /* count the cones in the union */
290 for (p = pStart; p < pEnd; p++) {
291 PetscInt dof, uOff;
292
293 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr);
294 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr);
295 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr);
296 coneSizes[uOff] = dof;
297 }
298 for (p = pRefStart; p < pRefEnd; p++) {
299 PetscInt dof, uDof, uOff;
300
301 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr);
302 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr);
303 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr);
304 if (uDof) {
305 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr);
306 coneSizes[uOff] = dof;
307 }
308 }
309 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr);
310 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr);
311 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr);
312 /* write the cones in the union */
313 for (p = pStart; p < pEnd; p++) {
314 PetscInt dof, uOff, c, cOff;
315 const PetscInt *cone, *orientation;
316
317 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr);
318 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr);
319 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr);
320 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr);
321 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr);
322 for (c = 0; c < dof; c++) {
323 PetscInt e, eOff;
324 e = cone[c];
325 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr);
326 unionCones[cOff + c] = eOff;
327 unionOrientations[cOff + c] = orientation[c];
328 }
329 }
330 for (p = pRefStart; p < pRefEnd; p++) {
331 PetscInt dof, uDof, uOff, c, cOff;
332 const PetscInt *cone, *orientation;
333
334 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr);
335 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr);
336 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr);
337 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr);
338 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr);
339 if (uDof) {
340 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr);
341 for (c = 0; c < dof; c++) {
342 PetscInt e, eOff, eDof;
343
344 e = cone[c];
345 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr);
346 if (eDof) {
347 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr);
348 }
349 else {
350 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr);
351 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr);
352 }
353 unionCones[cOff + c] = eOff;
354 unionOrientations[cOff + c] = orientation[c];
355 }
356 }
357 }
358 /* get the coordinates */
359 {
360 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
361 PetscSection KcoordsSec, KrefCoordsSec;
362 Vec KcoordsVec, KrefCoordsVec;
363 PetscScalar *Kcoords;
364
365 ierr = DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr);
366 ierr = DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr);
367 ierr = DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr);
368 ierr = DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr);
369
370 numVerts = numDimPoints[0];
371 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr);
372 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr);
373
374 offset = 0;
375 for (v = vStart; v < vEnd; v++) {
376 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr);
377 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr);
378 for (d = 0; d < dim; d++) {
379 unionCoords[offset * dim + d] = Kcoords[d];
380 }
381 offset++;
382 }
383 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr);
384 for (v = vRefStart; v < vRefEnd; v++) {
385 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr);
386 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr);
387 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr);
388 if (vDof) {
389 for (d = 0; d < dim; d++) {
390 unionCoords[offset * dim + d] = Kcoords[d];
391 }
392 offset++;
393 }
394 }
395 }
396 ierr = DMCreate(comm,ref);CHKERRQ(ierr);
397 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr);
398 ierr = DMSetDimension(*ref,dim);CHKERRQ(ierr);
399 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr);
400 /* set the tree */
401 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr);
402 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr);
403 for (p = pRefStart; p < pRefEnd; p++) {
404 PetscInt uDof, uOff;
405
406 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr);
407 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr);
408 if (uDof) {
409 ierr = PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr);
410 }
411 }
412 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr);
413 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr);
414 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr);
415 for (p = pRefStart; p < pRefEnd; p++) {
416 PetscInt uDof, uOff;
417
418 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr);
419 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr);
420 if (uDof) {
421 PetscInt pOff, parent, parentU;
422 ierr = PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr);
423 ierr = DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr);
424 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr);
425 parents[pOff] = parentU;
426 childIDs[pOff] = uOff;
427 }
428 }
429 ierr = DMPlexCreateReferenceTree_SetTree(*ref,parentSection,parents,childIDs);CHKERRQ(ierr);
430 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr);
431 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr);
432
433 /* clean up */
434 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr);
435 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr);
436 ierr = ISDestroy(&perm);CHKERRQ(ierr);
437 ierr = PetscFree(unionCoords);CHKERRQ(ierr);
438 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr);
439 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr);
440 PetscFunctionReturn(0);
441 }
442
443 /*@
444 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.
445
446 Collective
447
448 Input Parameters:
449 + comm - the MPI communicator
450 . dim - the spatial dimension
451 - simplex - Flag for simplex, otherwise use a tensor-product cell
452
453 Output Parameters:
454 . ref - the reference tree DMPlex object
455
456 Level: intermediate
457
458 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
459 @*/
DMPlexCreateDefaultReferenceTree(MPI_Comm comm,PetscInt dim,PetscBool simplex,DM * ref)460 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
461 {
462 DM_Plex *mesh;
463 DM K, Kref;
464 PetscInt p, pStart, pEnd;
465 DMLabel identity;
466 PetscErrorCode ierr;
467
468 PetscFunctionBegin;
469 #if 1
470 comm = PETSC_COMM_SELF;
471 #endif
472 /* create a reference element */
473 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr);
474 ierr = DMCreateLabel(K, "identity");CHKERRQ(ierr);
475 ierr = DMGetLabel(K, "identity", &identity);CHKERRQ(ierr);
476 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr);
477 for (p = pStart; p < pEnd; p++) {
478 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr);
479 }
480 /* refine it */
481 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr);
482
483 /* the reference tree is the union of these two, without duplicating
484 * points that appear in both */
485 ierr = DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);CHKERRQ(ierr);
486 mesh = (DM_Plex *) (*ref)->data;
487 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
488 ierr = DMDestroy(&K);CHKERRQ(ierr);
489 ierr = DMDestroy(&Kref);CHKERRQ(ierr);
490 PetscFunctionReturn(0);
491 }
492
DMPlexTreeSymmetrize(DM dm)493 static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
494 {
495 DM_Plex *mesh = (DM_Plex *)dm->data;
496 PetscSection childSec, pSec;
497 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
498 PetscInt *offsets, *children, pStart, pEnd;
499 PetscErrorCode ierr;
500
501 PetscFunctionBegin;
502 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
503 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr);
504 ierr = PetscFree(mesh->children);CHKERRQ(ierr);
505 pSec = mesh->parentSection;
506 if (!pSec) PetscFunctionReturn(0);
507 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr);
508 for (p = 0; p < pSize; p++) {
509 PetscInt par = mesh->parents[p];
510
511 parMax = PetscMax(parMax,par+1);
512 parMin = PetscMin(parMin,par);
513 }
514 if (parMin > parMax) {
515 parMin = -1;
516 parMax = -1;
517 }
518 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr);
519 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr);
520 for (p = 0; p < pSize; p++) {
521 PetscInt par = mesh->parents[p];
522
523 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr);
524 }
525 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr);
526 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr);
527 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr);
528 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr);
529 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr);
530 for (p = pStart; p < pEnd; p++) {
531 PetscInt dof, off, i;
532
533 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr);
534 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr);
535 for (i = 0; i < dof; i++) {
536 PetscInt par = mesh->parents[off + i], cOff;
537
538 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr);
539 children[cOff + offsets[par-parMin]++] = p;
540 }
541 }
542 mesh->childSection = childSec;
543 mesh->children = children;
544 ierr = PetscFree(offsets);CHKERRQ(ierr);
545 PetscFunctionReturn(0);
546 }
547
AnchorsFlatten(PetscSection section,IS is,PetscSection * sectionNew,IS * isNew)548 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
549 {
550 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
551 const PetscInt *vals;
552 PetscSection secNew;
553 PetscBool anyNew, globalAnyNew;
554 PetscBool compress;
555 PetscErrorCode ierr;
556
557 PetscFunctionBegin;
558 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr);
559 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr);
560 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr);
561 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr);
562 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr);
563 for (i = 0; i < size; i++) {
564 PetscInt dof;
565
566 p = vals[i];
567 if (p < pStart || p >= pEnd) continue;
568 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
569 if (dof) break;
570 }
571 if (i == size) {
572 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr);
573 anyNew = PETSC_FALSE;
574 compress = PETSC_FALSE;
575 sizeNew = 0;
576 }
577 else {
578 anyNew = PETSC_TRUE;
579 for (p = pStart; p < pEnd; p++) {
580 PetscInt dof, off;
581
582 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
583 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
584 for (i = 0; i < dof; i++) {
585 PetscInt q = vals[off + i], qDof = 0;
586
587 if (q >= pStart && q < pEnd) {
588 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr);
589 }
590 if (qDof) {
591 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr);
592 }
593 else {
594 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr);
595 }
596 }
597 }
598 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr);
599 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr);
600 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr);
601 compress = PETSC_FALSE;
602 for (p = pStart; p < pEnd; p++) {
603 PetscInt dof, off, count, offNew, dofNew;
604
605 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
606 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
607 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr);
608 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr);
609 count = 0;
610 for (i = 0; i < dof; i++) {
611 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
612
613 if (q >= pStart && q < pEnd) {
614 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr);
615 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr);
616 }
617 if (qDof) {
618 PetscInt oldCount = count;
619
620 for (j = 0; j < qDof; j++) {
621 PetscInt k, r = vals[qOff + j];
622
623 for (k = 0; k < oldCount; k++) {
624 if (valsNew[offNew + k] == r) {
625 break;
626 }
627 }
628 if (k == oldCount) {
629 valsNew[offNew + count++] = r;
630 }
631 }
632 }
633 else {
634 PetscInt k, oldCount = count;
635
636 for (k = 0; k < oldCount; k++) {
637 if (valsNew[offNew + k] == q) {
638 break;
639 }
640 }
641 if (k == oldCount) {
642 valsNew[offNew + count++] = q;
643 }
644 }
645 }
646 if (count < dofNew) {
647 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr);
648 compress = PETSC_TRUE;
649 }
650 }
651 }
652 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr);
653 ierr = MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr);
654 if (!globalAnyNew) {
655 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr);
656 *sectionNew = NULL;
657 *isNew = NULL;
658 }
659 else {
660 PetscBool globalCompress;
661
662 ierr = MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr);
663 if (compress) {
664 PetscSection secComp;
665 PetscInt *valsComp = NULL;
666
667 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr);
668 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr);
669 for (p = pStart; p < pEnd; p++) {
670 PetscInt dof;
671
672 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr);
673 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr);
674 }
675 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr);
676 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr);
677 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr);
678 for (p = pStart; p < pEnd; p++) {
679 PetscInt dof, off, offNew, j;
680
681 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr);
682 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr);
683 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr);
684 for (j = 0; j < dof; j++) {
685 valsComp[offNew + j] = valsNew[off + j];
686 }
687 }
688 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr);
689 secNew = secComp;
690 ierr = PetscFree(valsNew);CHKERRQ(ierr);
691 valsNew = valsComp;
692 }
693 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr);
694 }
695 PetscFunctionReturn(0);
696 }
697
DMPlexCreateAnchors_Tree(DM dm)698 static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
699 {
700 PetscInt p, pStart, pEnd, *anchors, size;
701 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
702 PetscSection aSec;
703 DMLabel canonLabel;
704 IS aIS;
705 PetscErrorCode ierr;
706
707 PetscFunctionBegin;
708 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
709 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
710 ierr = DMGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr);
711 for (p = pStart; p < pEnd; p++) {
712 PetscInt parent;
713
714 if (canonLabel) {
715 PetscInt canon;
716
717 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr);
718 if (p != canon) continue;
719 }
720 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr);
721 if (parent != p) {
722 aMin = PetscMin(aMin,p);
723 aMax = PetscMax(aMax,p+1);
724 }
725 }
726 if (aMin > aMax) {
727 aMin = -1;
728 aMax = -1;
729 }
730 ierr = PetscSectionCreate(PETSC_COMM_SELF,&aSec);CHKERRQ(ierr);
731 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr);
732 for (p = aMin; p < aMax; p++) {
733 PetscInt parent, ancestor = p;
734
735 if (canonLabel) {
736 PetscInt canon;
737
738 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr);
739 if (p != canon) continue;
740 }
741 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr);
742 while (parent != ancestor) {
743 ancestor = parent;
744 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr);
745 }
746 if (ancestor != p) {
747 PetscInt closureSize, *closure = NULL;
748
749 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
750 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr);
751 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
752 }
753 }
754 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr);
755 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr);
756 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr);
757 for (p = aMin; p < aMax; p++) {
758 PetscInt parent, ancestor = p;
759
760 if (canonLabel) {
761 PetscInt canon;
762
763 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr);
764 if (p != canon) continue;
765 }
766 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr);
767 while (parent != ancestor) {
768 ancestor = parent;
769 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr);
770 }
771 if (ancestor != p) {
772 PetscInt j, closureSize, *closure = NULL, aOff;
773
774 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr);
775
776 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
777 for (j = 0; j < closureSize; j++) {
778 anchors[aOff + j] = closure[2*j];
779 }
780 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
781 }
782 }
783 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr);
784 {
785 PetscSection aSecNew = aSec;
786 IS aISNew = aIS;
787
788 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr);
789 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr);
790 while (aSecNew) {
791 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr);
792 ierr = ISDestroy(&aIS);CHKERRQ(ierr);
793 aSec = aSecNew;
794 aIS = aISNew;
795 aSecNew = NULL;
796 aISNew = NULL;
797 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr);
798 }
799 }
800 ierr = DMPlexSetAnchors(dm,aSec,aIS);CHKERRQ(ierr);
801 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr);
802 ierr = ISDestroy(&aIS);CHKERRQ(ierr);
803 PetscFunctionReturn(0);
804 }
805
DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt * dof,PetscInt * numTrueSupp)806 static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
807 {
808 PetscErrorCode ierr;
809
810 PetscFunctionBegin;
811 if (numTrueSupp[p] == -1) {
812 PetscInt i, alldof;
813 const PetscInt *supp;
814 PetscInt count = 0;
815
816 ierr = DMPlexGetSupportSize(dm,p,&alldof);CHKERRQ(ierr);
817 ierr = DMPlexGetSupport(dm,p,&supp);CHKERRQ(ierr);
818 for (i = 0; i < alldof; i++) {
819 PetscInt q = supp[i], numCones, j;
820 const PetscInt *cone;
821
822 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr);
823 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr);
824 for (j = 0; j < numCones; j++) {
825 if (cone[j] == p) break;
826 }
827 if (j < numCones) count++;
828 }
829 numTrueSupp[p] = count;
830 }
831 *dof = numTrueSupp[p];
832 PetscFunctionReturn(0);
833 }
834
DMPlexTreeExchangeSupports(DM dm)835 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
836 {
837 DM_Plex *mesh = (DM_Plex *)dm->data;
838 PetscSection newSupportSection;
839 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
840 PetscInt *numTrueSupp;
841 PetscInt *offsets;
842 PetscErrorCode ierr;
843
844 PetscFunctionBegin;
845 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
846 /* symmetrize the hierarchy */
847 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr);
848 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);CHKERRQ(ierr);
849 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
850 ierr = PetscSectionSetChart(newSupportSection,pStart,pEnd);CHKERRQ(ierr);
851 ierr = PetscCalloc1(pEnd,&offsets);CHKERRQ(ierr);
852 ierr = PetscMalloc1(pEnd,&numTrueSupp);CHKERRQ(ierr);
853 for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
854 /* if a point is in the (true) support of q, it should be in the support of
855 * parent(q) */
856 for (d = 0; d <= depth; d++) {
857 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr);
858 for (p = pStart; p < pEnd; ++p) {
859 PetscInt dof, q, qdof, parent;
860
861 ierr = DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);CHKERRQ(ierr);
862 ierr = PetscSectionAddDof(newSupportSection, p, dof);CHKERRQ(ierr);
863 q = p;
864 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr);
865 while (parent != q && parent >= pStart && parent < pEnd) {
866 q = parent;
867
868 ierr = DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);CHKERRQ(ierr);
869 ierr = PetscSectionAddDof(newSupportSection,p,qdof);CHKERRQ(ierr);
870 ierr = PetscSectionAddDof(newSupportSection,q,dof);CHKERRQ(ierr);
871 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr);
872 }
873 }
874 }
875 ierr = PetscSectionSetUp(newSupportSection);CHKERRQ(ierr);
876 ierr = PetscSectionGetStorageSize(newSupportSection,&newSize);CHKERRQ(ierr);
877 ierr = PetscMalloc1(newSize,&newSupports);CHKERRQ(ierr);
878 for (d = 0; d <= depth; d++) {
879 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr);
880 for (p = pStart; p < pEnd; p++) {
881 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
882
883 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr);
884 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr);
885 ierr = PetscSectionGetDof(newSupportSection, p, &newDof);CHKERRQ(ierr);
886 ierr = PetscSectionGetOffset(newSupportSection, p, &newOff);CHKERRQ(ierr);
887 for (i = 0; i < dof; i++) {
888 PetscInt numCones, j;
889 const PetscInt *cone;
890 PetscInt q = mesh->supports[off + i];
891
892 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr);
893 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr);
894 for (j = 0; j < numCones; j++) {
895 if (cone[j] == p) break;
896 }
897 if (j < numCones) newSupports[newOff+offsets[p]++] = q;
898 }
899 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);
900
901 q = p;
902 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr);
903 while (parent != q && parent >= pStart && parent < pEnd) {
904 q = parent;
905 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr);
906 ierr = PetscSectionGetOffset(mesh->supportSection, q, &qoff);CHKERRQ(ierr);
907 ierr = PetscSectionGetOffset(newSupportSection, q, &newqOff);CHKERRQ(ierr);
908 for (i = 0; i < qdof; i++) {
909 PetscInt numCones, j;
910 const PetscInt *cone;
911 PetscInt r = mesh->supports[qoff + i];
912
913 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr);
914 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr);
915 for (j = 0; j < numCones; j++) {
916 if (cone[j] == q) break;
917 }
918 if (j < numCones) newSupports[newOff+offsets[p]++] = r;
919 }
920 for (i = 0; i < dof; i++) {
921 PetscInt numCones, j;
922 const PetscInt *cone;
923 PetscInt r = mesh->supports[off + i];
924
925 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr);
926 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr);
927 for (j = 0; j < numCones; j++) {
928 if (cone[j] == p) break;
929 }
930 if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
931 }
932 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr);
933 }
934 }
935 }
936 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr);
937 mesh->supportSection = newSupportSection;
938 ierr = PetscFree(mesh->supports);CHKERRQ(ierr);
939 mesh->supports = newSupports;
940 ierr = PetscFree(offsets);CHKERRQ(ierr);
941 ierr = PetscFree(numTrueSupp);CHKERRQ(ierr);
942
943 PetscFunctionReturn(0);
944 }
945
946 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
947 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);
948
DMPlexSetTree_Internal(DM dm,PetscSection parentSection,PetscInt * parents,PetscInt * childIDs,PetscBool computeCanonical,PetscBool exchangeSupports)949 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
950 {
951 DM_Plex *mesh = (DM_Plex *)dm->data;
952 DM refTree;
953 PetscInt size;
954 PetscErrorCode ierr;
955
956 PetscFunctionBegin;
957 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
958 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2);
959 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr);
960 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr);
961 mesh->parentSection = parentSection;
962 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr);
963 if (parents != mesh->parents) {
964 ierr = PetscFree(mesh->parents);CHKERRQ(ierr);
965 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr);
966 ierr = PetscArraycpy(mesh->parents, parents, size);CHKERRQ(ierr);
967 }
968 if (childIDs != mesh->childIDs) {
969 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr);
970 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr);
971 ierr = PetscArraycpy(mesh->childIDs, childIDs, size);CHKERRQ(ierr);
972 }
973 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr);
974 if (refTree) {
975 DMLabel canonLabel;
976
977 ierr = DMGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr);
978 if (canonLabel) {
979 PetscInt i;
980
981 for (i = 0; i < size; i++) {
982 PetscInt canon;
983 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr);
984 if (canon >= 0) {
985 mesh->childIDs[i] = canon;
986 }
987 }
988 }
989 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
990 } else {
991 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
992 }
993 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr);
994 if (computeCanonical) {
995 PetscInt d, dim;
996
997 /* add the canonical label */
998 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
999 ierr = DMCreateLabel(dm,"canonical");CHKERRQ(ierr);
1000 for (d = 0; d <= dim; d++) {
1001 PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
1002 const PetscInt *cChildren;
1003
1004 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr);
1005 for (p = dStart; p < dEnd; p++) {
1006 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr);
1007 if (cNumChildren) {
1008 canon = p;
1009 break;
1010 }
1011 }
1012 if (canon == -1) continue;
1013 for (p = dStart; p < dEnd; p++) {
1014 PetscInt numChildren, i;
1015 const PetscInt *children;
1016
1017 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr);
1018 if (numChildren) {
1019 if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren);
1020 ierr = DMSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr);
1021 for (i = 0; i < numChildren; i++) {
1022 ierr = DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr);
1023 }
1024 }
1025 }
1026 }
1027 }
1028 if (exchangeSupports) {
1029 ierr = DMPlexTreeExchangeSupports(dm);CHKERRQ(ierr);
1030 }
1031 mesh->createanchors = DMPlexCreateAnchors_Tree;
1032 /* reset anchors */
1033 ierr = DMPlexSetAnchors(dm,NULL,NULL);CHKERRQ(ierr);
1034 PetscFunctionReturn(0);
1035 }
1036
1037 /*@
1038 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
1039 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its
1040 tree root.
1041
1042 Collective on dm
1043
1044 Input Parameters:
1045 + dm - the DMPlex object
1046 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1047 offset indexes the parent and childID list; the reference count of parentSection is incremented
1048 . parents - a list of the point parents; copied, can be destroyed
1049 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1050 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
1051
1052 Level: intermediate
1053
1054 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1055 @*/
DMPlexSetTree(DM dm,PetscSection parentSection,PetscInt parents[],PetscInt childIDs[])1056 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1057 {
1058 PetscErrorCode ierr;
1059
1060 PetscFunctionBegin;
1061 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr);
1062 PetscFunctionReturn(0);
1063 }
1064
1065 /*@
1066 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1067 Collective on dm
1068
1069 Input Parameters:
1070 . dm - the DMPlex object
1071
1072 Output Parameters:
1073 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1074 offset indexes the parent and childID list
1075 . parents - a list of the point parents
1076 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1077 the child corresponds to the point in the reference tree with index childID
1078 . childSection - the inverse of the parent section
1079 - children - a list of the point children
1080
1081 Level: intermediate
1082
1083 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1084 @*/
DMPlexGetTree(DM dm,PetscSection * parentSection,PetscInt * parents[],PetscInt * childIDs[],PetscSection * childSection,PetscInt * children[])1085 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1086 {
1087 DM_Plex *mesh = (DM_Plex *)dm->data;
1088
1089 PetscFunctionBegin;
1090 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1091 if (parentSection) *parentSection = mesh->parentSection;
1092 if (parents) *parents = mesh->parents;
1093 if (childIDs) *childIDs = mesh->childIDs;
1094 if (childSection) *childSection = mesh->childSection;
1095 if (children) *children = mesh->children;
1096 PetscFunctionReturn(0);
1097 }
1098
1099 /*@
1100 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG)
1101
1102 Input Parameters:
1103 + dm - the DMPlex object
1104 - point - the query point
1105
1106 Output Parameters:
1107 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1108 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1109 does not have a parent
1110
1111 Level: intermediate
1112
1113 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1114 @*/
DMPlexGetTreeParent(DM dm,PetscInt point,PetscInt * parent,PetscInt * childID)1115 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1116 {
1117 DM_Plex *mesh = (DM_Plex *)dm->data;
1118 PetscSection pSec;
1119 PetscErrorCode ierr;
1120
1121 PetscFunctionBegin;
1122 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1123 pSec = mesh->parentSection;
1124 if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1125 PetscInt dof;
1126
1127 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr);
1128 if (dof) {
1129 PetscInt off;
1130
1131 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr);
1132 if (parent) *parent = mesh->parents[off];
1133 if (childID) *childID = mesh->childIDs[off];
1134 PetscFunctionReturn(0);
1135 }
1136 }
1137 if (parent) {
1138 *parent = point;
1139 }
1140 if (childID) {
1141 *childID = 0;
1142 }
1143 PetscFunctionReturn(0);
1144 }
1145
1146 /*@C
1147 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG)
1148
1149 Input Parameters:
1150 + dm - the DMPlex object
1151 - point - the query point
1152
1153 Output Parameters:
1154 + numChildren - if not NULL, set to the number of children
1155 - children - if not NULL, set to a list children, or set to NULL if the point has no children
1156
1157 Level: intermediate
1158
1159 Fortran Notes:
1160 Since it returns an array, this routine is only available in Fortran 90, and you must
1161 include petsc.h90 in your code.
1162
1163 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1164 @*/
DMPlexGetTreeChildren(DM dm,PetscInt point,PetscInt * numChildren,const PetscInt * children[])1165 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1166 {
1167 DM_Plex *mesh = (DM_Plex *)dm->data;
1168 PetscSection childSec;
1169 PetscInt dof = 0;
1170 PetscErrorCode ierr;
1171
1172 PetscFunctionBegin;
1173 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1174 childSec = mesh->childSection;
1175 if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1176 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr);
1177 }
1178 if (numChildren) *numChildren = dof;
1179 if (children) {
1180 if (dof) {
1181 PetscInt off;
1182
1183 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr);
1184 *children = &mesh->children[off];
1185 }
1186 else {
1187 *children = NULL;
1188 }
1189 }
1190 PetscFunctionReturn(0);
1191 }
1192
EvaluateBasis(PetscSpace space,PetscInt nBasis,PetscInt nFunctionals,PetscInt nComps,PetscInt nPoints,const PetscInt * pointsPerFn,const PetscReal * points,const PetscReal * weights,PetscReal * work,Mat basisAtPoints)1193 static PetscErrorCode EvaluateBasis(PetscSpace space, PetscInt nBasis, PetscInt nFunctionals, PetscInt nComps, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints)
1194 {
1195 PetscInt f, b, p, c, offset, qPoints;
1196 PetscErrorCode ierr;
1197
1198 PetscFunctionBegin;
1199 ierr = PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);CHKERRQ(ierr);
1200 for (f = 0, offset = 0; f < nFunctionals; f++) {
1201 qPoints = pointsPerFn[f];
1202 for (b = 0; b < nBasis; b++) {
1203 PetscScalar val = 0.;
1204
1205 for (p = 0; p < qPoints; p++) {
1206 for (c = 0; c < nComps; c++) {
1207 val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1208 }
1209 }
1210 ierr = MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES);CHKERRQ(ierr);
1211 }
1212 offset += qPoints;
1213 }
1214 ierr = MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1215 ierr = MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1216 PetscFunctionReturn(0);
1217 }
1218
DMPlexComputeAnchorMatrix_Tree_Direct(DM dm,PetscSection section,PetscSection cSec,Mat cMat)1219 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1220 {
1221 PetscDS ds;
1222 PetscInt spdim;
1223 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1224 const PetscInt *anchors;
1225 PetscSection aSec;
1226 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1227 IS aIS;
1228 PetscErrorCode ierr;
1229
1230 PetscFunctionBegin;
1231 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
1232 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr);
1233 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
1234 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
1235 ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr);
1236 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr);
1237 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr);
1238 ierr = DMGetDimension(dm,&spdim);CHKERRQ(ierr);
1239 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr);
1240
1241 for (f = 0; f < numFields; f++) {
1242 PetscObject disc;
1243 PetscClassId id;
1244 PetscSpace bspace;
1245 PetscDualSpace dspace;
1246 PetscInt i, j, k, nPoints, Nc, offset;
1247 PetscInt fSize, maxDof;
1248 PetscReal *weights, *pointsRef, *pointsReal, *work;
1249 PetscScalar *scwork;
1250 const PetscScalar *X;
1251 PetscInt *sizes, *workIndRow, *workIndCol;
1252 Mat Amat, Bmat, Xmat;
1253 const PetscInt *numDof = NULL;
1254 const PetscInt ***perms = NULL;
1255 const PetscScalar ***flips = NULL;
1256
1257 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr);
1258 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr);
1259 if (id == PETSCFE_CLASSID) {
1260 PetscFE fe = (PetscFE) disc;
1261
1262 ierr = PetscFEGetBasisSpace(fe,&bspace);CHKERRQ(ierr);
1263 ierr = PetscFEGetDualSpace(fe,&dspace);CHKERRQ(ierr);
1264 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr);
1265 ierr = PetscFEGetNumComponents(fe,&Nc);CHKERRQ(ierr);
1266 }
1267 else if (id == PETSCFV_CLASSID) {
1268 PetscFV fv = (PetscFV) disc;
1269
1270 ierr = PetscFVGetNumComponents(fv,&Nc);CHKERRQ(ierr);
1271 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);CHKERRQ(ierr);
1272 ierr = PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr);
1273 ierr = PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE);CHKERRQ(ierr);
1274 ierr = PetscSpaceSetNumComponents(bspace,Nc);CHKERRQ(ierr);
1275 ierr = PetscSpaceSetNumVariables(bspace,spdim);CHKERRQ(ierr);
1276 ierr = PetscSpaceSetUp(bspace);CHKERRQ(ierr);
1277 ierr = PetscFVGetDualSpace(fv,&dspace);CHKERRQ(ierr);
1278 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr);
1279 }
1280 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1281 ierr = PetscDualSpaceGetNumDof(dspace,&numDof);CHKERRQ(ierr);
1282 for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);}
1283 ierr = PetscDualSpaceGetSymmetries(dspace,&perms,&flips);CHKERRQ(ierr);
1284
1285 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr);
1286 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr);
1287 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr);
1288 ierr = MatSetUp(Amat);CHKERRQ(ierr);
1289 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr);
1290 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr);
1291 nPoints = 0;
1292 for (i = 0; i < fSize; i++) {
1293 PetscInt qPoints, thisNc;
1294 PetscQuadrature quad;
1295
1296 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr);
1297 ierr = PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL);CHKERRQ(ierr);
1298 if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
1299 nPoints += qPoints;
1300 }
1301 ierr = PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol);CHKERRQ(ierr);
1302 ierr = PetscMalloc1(maxDof * maxDof,&scwork);CHKERRQ(ierr);
1303 offset = 0;
1304 for (i = 0; i < fSize; i++) {
1305 PetscInt qPoints;
1306 const PetscReal *p, *w;
1307 PetscQuadrature quad;
1308
1309 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr);
1310 ierr = PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w);CHKERRQ(ierr);
1311 ierr = PetscArraycpy(weights+Nc*offset,w,Nc*qPoints);CHKERRQ(ierr);
1312 ierr = PetscArraycpy(pointsRef+spdim*offset,p,spdim*qPoints);CHKERRQ(ierr);
1313 sizes[i] = qPoints;
1314 offset += qPoints;
1315 }
1316 ierr = EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat);CHKERRQ(ierr);
1317 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr);
1318 for (c = cStart; c < cEnd; c++) {
1319 PetscInt parent;
1320 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1321 PetscInt *childOffsets, *parentOffsets;
1322
1323 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr);
1324 if (parent == c) continue;
1325 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1326 for (i = 0; i < closureSize; i++) {
1327 PetscInt p = closure[2*i];
1328 PetscInt conDof;
1329
1330 if (p < conStart || p >= conEnd) continue;
1331 if (numFields) {
1332 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr);
1333 }
1334 else {
1335 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr);
1336 }
1337 if (conDof) break;
1338 }
1339 if (i == closureSize) {
1340 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1341 continue;
1342 }
1343
1344 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr);
1345 ierr = DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr);
1346 for (i = 0; i < nPoints; i++) {
1347 const PetscReal xi0[3] = {-1.,-1.,-1.};
1348
1349 CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp);
1350 CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1351 }
1352 ierr = EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);CHKERRQ(ierr);
1353 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr);
1354 ierr = MatDenseGetArrayRead(Xmat,&X);CHKERRQ(ierr);
1355 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr);
1356 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr);
1357 childOffsets[0] = 0;
1358 for (i = 0; i < closureSize; i++) {
1359 PetscInt p = closure[2*i];
1360 PetscInt dof;
1361
1362 if (numFields) {
1363 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr);
1364 }
1365 else {
1366 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr);
1367 }
1368 childOffsets[i+1]=childOffsets[i]+dof;
1369 }
1370 parentOffsets[0] = 0;
1371 for (i = 0; i < closureSizeP; i++) {
1372 PetscInt p = closureP[2*i];
1373 PetscInt dof;
1374
1375 if (numFields) {
1376 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr);
1377 }
1378 else {
1379 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr);
1380 }
1381 parentOffsets[i+1]=parentOffsets[i]+dof;
1382 }
1383 for (i = 0; i < closureSize; i++) {
1384 PetscInt conDof, conOff, aDof, aOff, nWork;
1385 PetscInt p = closure[2*i];
1386 PetscInt o = closure[2*i+1];
1387 const PetscInt *perm;
1388 const PetscScalar *flip;
1389
1390 if (p < conStart || p >= conEnd) continue;
1391 if (numFields) {
1392 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr);
1393 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr);
1394 }
1395 else {
1396 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr);
1397 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr);
1398 }
1399 if (!conDof) continue;
1400 perm = (perms && perms[i]) ? perms[i][o] : NULL;
1401 flip = (flips && flips[i]) ? flips[i][o] : NULL;
1402 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr);
1403 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr);
1404 nWork = childOffsets[i+1]-childOffsets[i];
1405 for (k = 0; k < aDof; k++) {
1406 PetscInt a = anchors[aOff + k];
1407 PetscInt aSecDof, aSecOff;
1408
1409 if (numFields) {
1410 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr);
1411 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr);
1412 }
1413 else {
1414 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr);
1415 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr);
1416 }
1417 if (!aSecDof) continue;
1418
1419 for (j = 0; j < closureSizeP; j++) {
1420 PetscInt q = closureP[2*j];
1421 PetscInt oq = closureP[2*j+1];
1422
1423 if (q == a) {
1424 PetscInt r, s, nWorkP;
1425 const PetscInt *permP;
1426 const PetscScalar *flipP;
1427
1428 permP = (perms && perms[j]) ? perms[j][oq] : NULL;
1429 flipP = (flips && flips[j]) ? flips[j][oq] : NULL;
1430 nWorkP = parentOffsets[j+1]-parentOffsets[j];
1431 /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1432 * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
1433 * column-major, so transpose-transpose = do nothing */
1434 for (r = 0; r < nWork; r++) {
1435 for (s = 0; s < nWorkP; s++) {
1436 scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1437 }
1438 }
1439 for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + r;}
1440 for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;}
1441 if (flip) {
1442 for (r = 0; r < nWork; r++) {
1443 for (s = 0; s < nWorkP; s++) {
1444 scwork[r * nWorkP + s] *= flip[r];
1445 }
1446 }
1447 }
1448 if (flipP) {
1449 for (r = 0; r < nWork; r++) {
1450 for (s = 0; s < nWorkP; s++) {
1451 scwork[r * nWorkP + s] *= flipP[s];
1452 }
1453 }
1454 }
1455 ierr = MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);CHKERRQ(ierr);
1456 break;
1457 }
1458 }
1459 }
1460 }
1461 ierr = MatDenseRestoreArrayRead(Xmat,&X);CHKERRQ(ierr);
1462 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr);
1463 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1464 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr);
1465 }
1466 ierr = MatDestroy(&Amat);CHKERRQ(ierr);
1467 ierr = MatDestroy(&Bmat);CHKERRQ(ierr);
1468 ierr = MatDestroy(&Xmat);CHKERRQ(ierr);
1469 ierr = PetscFree(scwork);CHKERRQ(ierr);
1470 ierr = PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);CHKERRQ(ierr);
1471 if (id == PETSCFV_CLASSID) {
1472 ierr = PetscSpaceDestroy(&bspace);CHKERRQ(ierr);
1473 }
1474 }
1475 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1476 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1477 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr);
1478 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr);
1479
1480 PetscFunctionReturn(0);
1481 }
1482
DMPlexReferenceTreeGetChildrenMatrices(DM refTree,PetscScalar **** childrenMats,PetscInt *** childrenN)1483 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1484 {
1485 Mat refCmat;
1486 PetscDS ds;
1487 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1488 PetscScalar ***refPointFieldMats;
1489 PetscSection refConSec, refAnSec, refSection;
1490 IS refAnIS;
1491 const PetscInt *refAnchors;
1492 const PetscInt **perms;
1493 const PetscScalar **flips;
1494 PetscErrorCode ierr;
1495
1496 PetscFunctionBegin;
1497 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr);
1498 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
1499 maxFields = PetscMax(1,numFields);
1500 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr);
1501 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr);
1502 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr);
1503 ierr = DMGetLocalSection(refTree,&refSection);CHKERRQ(ierr);
1504 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
1505 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr);
1506 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr);
1507 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr);
1508 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr);
1509 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr);
1510 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr);
1511 for (p = pRefStart; p < pRefEnd; p++) {
1512 PetscInt parent, closureSize, *closure = NULL, pDof;
1513
1514 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr);
1515 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr);
1516 if (!pDof || parent == p) continue;
1517
1518 ierr = PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr);
1519 ierr = PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr);
1520 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1521 for (f = 0; f < maxFields; f++) {
1522 PetscInt cDof, cOff, numCols, r, i;
1523
1524 if (f < numFields) {
1525 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr);
1526 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr);
1527 ierr = PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1528 } else {
1529 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr);
1530 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr);
1531 ierr = PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1532 }
1533
1534 for (r = 0; r < cDof; r++) {
1535 rows[r] = cOff + r;
1536 }
1537 numCols = 0;
1538 for (i = 0; i < closureSize; i++) {
1539 PetscInt q = closure[2*i];
1540 PetscInt aDof, aOff, j;
1541 const PetscInt *perm = perms ? perms[i] : NULL;
1542
1543 if (numFields) {
1544 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr);
1545 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr);
1546 }
1547 else {
1548 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr);
1549 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr);
1550 }
1551
1552 for (j = 0; j < aDof; j++) {
1553 cols[numCols++] = aOff + (perm ? perm[j] : j);
1554 }
1555 }
1556 refPointFieldN[p-pRefStart][f] = numCols;
1557 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr);
1558 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr);
1559 if (flips) {
1560 PetscInt colOff = 0;
1561
1562 for (i = 0; i < closureSize; i++) {
1563 PetscInt q = closure[2*i];
1564 PetscInt aDof, aOff, j;
1565 const PetscScalar *flip = flips ? flips[i] : NULL;
1566
1567 if (numFields) {
1568 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr);
1569 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr);
1570 }
1571 else {
1572 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr);
1573 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr);
1574 }
1575 if (flip) {
1576 PetscInt k;
1577 for (k = 0; k < cDof; k++) {
1578 for (j = 0; j < aDof; j++) {
1579 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j];
1580 }
1581 }
1582 }
1583 colOff += aDof;
1584 }
1585 }
1586 if (numFields) {
1587 ierr = PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1588 } else {
1589 ierr = PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1590 }
1591 }
1592 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1593 }
1594 *childrenMats = refPointFieldMats;
1595 *childrenN = refPointFieldN;
1596 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr);
1597 ierr = PetscFree(rows);CHKERRQ(ierr);
1598 ierr = PetscFree(cols);CHKERRQ(ierr);
1599 PetscFunctionReturn(0);
1600 }
1601
DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree,PetscScalar **** childrenMats,PetscInt *** childrenN)1602 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1603 {
1604 PetscDS ds;
1605 PetscInt **refPointFieldN;
1606 PetscScalar ***refPointFieldMats;
1607 PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f;
1608 PetscSection refConSec;
1609 PetscErrorCode ierr;
1610
1611 PetscFunctionBegin;
1612 refPointFieldN = *childrenN;
1613 *childrenN = NULL;
1614 refPointFieldMats = *childrenMats;
1615 *childrenMats = NULL;
1616 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr);
1617 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
1618 maxFields = PetscMax(1,numFields);
1619 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr);
1620 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
1621 for (p = pRefStart; p < pRefEnd; p++) {
1622 PetscInt parent, pDof;
1623
1624 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr);
1625 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr);
1626 if (!pDof || parent == p) continue;
1627
1628 for (f = 0; f < maxFields; f++) {
1629 PetscInt cDof;
1630
1631 if (numFields) {
1632 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr);
1633 }
1634 else {
1635 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr);
1636 }
1637
1638 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr);
1639 }
1640 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr);
1641 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr);
1642 }
1643 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr);
1644 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr);
1645 PetscFunctionReturn(0);
1646 }
1647
DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm,PetscSection section,PetscSection conSec,Mat cMat)1648 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1649 {
1650 DM refTree;
1651 PetscDS ds;
1652 Mat refCmat;
1653 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1654 PetscScalar ***refPointFieldMats, *pointWork;
1655 PetscSection refConSec, refAnSec, anSec;
1656 IS refAnIS, anIS;
1657 const PetscInt *anchors;
1658 PetscErrorCode ierr;
1659
1660 PetscFunctionBegin;
1661 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1662 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr);
1663 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
1664 maxFields = PetscMax(1,numFields);
1665 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr);
1666 ierr = DMCopyDisc(dm,refTree);CHKERRQ(ierr);
1667 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr);
1668 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr);
1669 ierr = DMPlexGetAnchors(dm,&anSec,&anIS);CHKERRQ(ierr);
1670 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr);
1671 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
1672 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr);
1673 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr);
1674 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr);
1675 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr);
1676
1677 /* step 1: get submats for every constrained point in the reference tree */
1678 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
1679
1680 /* step 2: compute the preorder */
1681 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
1682 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr);
1683 for (p = pStart; p < pEnd; p++) {
1684 perm[p - pStart] = p;
1685 iperm[p - pStart] = p-pStart;
1686 }
1687 for (p = 0; p < pEnd - pStart;) {
1688 PetscInt point = perm[p];
1689 PetscInt parent;
1690
1691 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr);
1692 if (parent == point) {
1693 p++;
1694 }
1695 else {
1696 PetscInt size, closureSize, *closure = NULL, i;
1697
1698 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1699 for (i = 0; i < closureSize; i++) {
1700 PetscInt q = closure[2*i];
1701 if (iperm[q-pStart] > iperm[point-pStart]) {
1702 /* swap */
1703 perm[p] = q;
1704 perm[iperm[q-pStart]] = point;
1705 iperm[point-pStart] = iperm[q-pStart];
1706 iperm[q-pStart] = p;
1707 break;
1708 }
1709 }
1710 size = closureSize;
1711 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1712 if (i == size) {
1713 p++;
1714 }
1715 }
1716 }
1717
1718 /* step 3: fill the constraint matrix */
1719 /* we are going to use a preorder progressive fill strategy. Mat doesn't
1720 * allow progressive fill without assembly, so we are going to set up the
1721 * values outside of the Mat first.
1722 */
1723 {
1724 PetscInt nRows, row, nnz;
1725 PetscBool done;
1726 const PetscInt *ia, *ja;
1727 PetscScalar *vals;
1728
1729 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr);
1730 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix");
1731 nnz = ia[nRows];
1732 /* malloc and then zero rows right before we fill them: this way valgrind
1733 * can tell if we are doing progressive fill in the wrong order */
1734 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr);
1735 for (p = 0; p < pEnd - pStart; p++) {
1736 PetscInt parent, childid, closureSize, *closure = NULL;
1737 PetscInt point = perm[p], pointDof;
1738
1739 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr);
1740 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1741 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr);
1742 if (!pointDof) continue;
1743 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1744 for (f = 0; f < maxFields; f++) {
1745 PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1746 PetscScalar *pointMat;
1747 const PetscInt **perms;
1748 const PetscScalar **flips;
1749
1750 if (numFields) {
1751 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr);
1752 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr);
1753 }
1754 else {
1755 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr);
1756 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr);
1757 }
1758 if (!cDof) continue;
1759 if (numFields) {ierr = PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);}
1760 else {ierr = PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr);}
1761
1762 /* make sure that every row for this point is the same size */
1763 if (PetscDefined(USE_DEBUG)) {
1764 for (r = 0; r < cDof; r++) {
1765 if (cDof > 1 && r) {
1766 if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %D vs. %D", (ia[cOff+r+1]-ia[cOff+r]), (ia[cOff+r]-ia[cOff+r-1]));
1767 }
1768 }
1769 }
1770 /* zero rows */
1771 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1772 vals[i] = 0.;
1773 }
1774 matOffset = ia[cOff];
1775 numFillCols = ia[cOff+1] - matOffset;
1776 pointMat = refPointFieldMats[childid-pRefStart][f];
1777 numCols = refPointFieldN[childid-pRefStart][f];
1778 offset = 0;
1779 for (i = 0; i < closureSize; i++) {
1780 PetscInt q = closure[2*i];
1781 PetscInt aDof, aOff, j, k, qConDof, qConOff;
1782 const PetscInt *perm = perms ? perms[i] : NULL;
1783 const PetscScalar *flip = flips ? flips[i] : NULL;
1784
1785 qConDof = qConOff = 0;
1786 if (numFields) {
1787 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr);
1788 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr);
1789 if (q >= conStart && q < conEnd) {
1790 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr);
1791 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr);
1792 }
1793 }
1794 else {
1795 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr);
1796 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr);
1797 if (q >= conStart && q < conEnd) {
1798 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr);
1799 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr);
1800 }
1801 }
1802 if (!aDof) continue;
1803 if (qConDof) {
1804 /* this point has anchors: its rows of the matrix should already
1805 * be filled, thanks to preordering */
1806 /* first multiply into pointWork, then set in matrix */
1807 PetscInt aMatOffset = ia[qConOff];
1808 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1809 for (r = 0; r < cDof; r++) {
1810 for (j = 0; j < aNumFillCols; j++) {
1811 PetscScalar inVal = 0;
1812 for (k = 0; k < aDof; k++) {
1813 PetscInt col = perm ? perm[k] : k;
1814
1815 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1816 }
1817 pointWork[r * aNumFillCols + j] = inVal;
1818 }
1819 }
1820 /* assume that the columns are sorted, spend less time searching */
1821 for (j = 0, k = 0; j < aNumFillCols; j++) {
1822 PetscInt col = ja[aMatOffset + j];
1823 for (;k < numFillCols; k++) {
1824 if (ja[matOffset + k] == col) {
1825 break;
1826 }
1827 }
1828 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1829 for (r = 0; r < cDof; r++) {
1830 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1831 }
1832 }
1833 }
1834 else {
1835 /* find where to put this portion of pointMat into the matrix */
1836 for (k = 0; k < numFillCols; k++) {
1837 if (ja[matOffset + k] == aOff) {
1838 break;
1839 }
1840 }
1841 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1842 for (r = 0; r < cDof; r++) {
1843 for (j = 0; j < aDof; j++) {
1844 PetscInt col = perm ? perm[j] : j;
1845
1846 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1847 }
1848 }
1849 }
1850 offset += aDof;
1851 }
1852 if (numFields) {
1853 ierr = PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1854 } else {
1855 ierr = PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr);
1856 }
1857 }
1858 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1859 }
1860 for (row = 0; row < nRows; row++) {
1861 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr);
1862 }
1863 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr);
1864 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1865 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1866 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1867 ierr = PetscFree(vals);CHKERRQ(ierr);
1868 }
1869
1870 /* clean up */
1871 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr);
1872 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr);
1873 ierr = PetscFree(pointWork);CHKERRQ(ierr);
1874 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
1875 PetscFunctionReturn(0);
1876 }
1877
1878 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1879 * a non-conforming mesh. Local refinement comes later */
DMPlexTreeRefineCell(DM dm,PetscInt cell,DM * ncdm)1880 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1881 {
1882 DM K;
1883 PetscMPIInt rank;
1884 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1885 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1886 PetscInt *Kembedding;
1887 PetscInt *cellClosure=NULL, nc;
1888 PetscScalar *newVertexCoords;
1889 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1890 PetscSection parentSection;
1891 PetscErrorCode ierr;
1892
1893 PetscFunctionBegin;
1894 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
1895 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1896 ierr = DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);CHKERRQ(ierr);
1897 ierr = DMSetDimension(*ncdm,dim);CHKERRQ(ierr);
1898
1899 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
1900 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);CHKERRQ(ierr);
1901 ierr = DMPlexGetReferenceTree(dm,&K);CHKERRQ(ierr);
1902 if (!rank) {
1903 /* compute the new charts */
1904 ierr = PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);CHKERRQ(ierr);
1905 offset = 0;
1906 for (d = 0; d <= dim; d++) {
1907 PetscInt pOldCount, kStart, kEnd, k;
1908
1909 pNewStart[d] = offset;
1910 ierr = DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);CHKERRQ(ierr);
1911 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr);
1912 pOldCount = pOldEnd[d] - pOldStart[d];
1913 /* adding the new points */
1914 pNewCount[d] = pOldCount + kEnd - kStart;
1915 if (!d) {
1916 /* removing the cell */
1917 pNewCount[d]--;
1918 }
1919 for (k = kStart; k < kEnd; k++) {
1920 PetscInt parent;
1921 ierr = DMPlexGetTreeParent(K,k,&parent,NULL);CHKERRQ(ierr);
1922 if (parent == k) {
1923 /* avoid double counting points that won't actually be new */
1924 pNewCount[d]--;
1925 }
1926 }
1927 pNewEnd[d] = pNewStart[d] + pNewCount[d];
1928 offset = pNewEnd[d];
1929
1930 }
1931 if (cell < pOldStart[0] || cell >= pOldEnd[0]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%d not in cell range [%d, %d)", cell, pOldStart[0], pOldEnd[0]);
1932 /* get the current closure of the cell that we are removing */
1933 ierr = DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr);
1934
1935 ierr = PetscMalloc1(pNewEnd[dim],&newConeSizes);CHKERRQ(ierr);
1936 {
1937 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1938
1939 ierr = DMPlexGetChart(K,&kStart,&kEnd);CHKERRQ(ierr);
1940 ierr = PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);CHKERRQ(ierr);
1941
1942 for (k = kStart; k < kEnd; k++) {
1943 perm[k - kStart] = k;
1944 iperm [k - kStart] = k - kStart;
1945 preOrient[k - kStart] = 0;
1946 }
1947
1948 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr);
1949 for (j = 1; j < closureSizeK; j++) {
1950 PetscInt parentOrientA = closureK[2*j+1];
1951 PetscInt parentOrientB = cellClosure[2*j+1];
1952 PetscInt p, q;
1953
1954 p = closureK[2*j];
1955 q = cellClosure[2*j];
1956 for (d = 0; d <= dim; d++) {
1957 if (q >= pOldStart[d] && q < pOldEnd[d]) {
1958 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1959 }
1960 }
1961 if (parentOrientA != parentOrientB) {
1962 PetscInt numChildren, i;
1963 const PetscInt *children;
1964
1965 ierr = DMPlexGetTreeChildren(K,p,&numChildren,&children);CHKERRQ(ierr);
1966 for (i = 0; i < numChildren; i++) {
1967 PetscInt kPerm, oPerm;
1968
1969 k = children[i];
1970 ierr = DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);CHKERRQ(ierr);
1971 /* perm = what refTree position I'm in */
1972 perm[kPerm-kStart] = k;
1973 /* iperm = who is at this position */
1974 iperm[k-kStart] = kPerm-kStart;
1975 preOrient[kPerm-kStart] = oPerm;
1976 }
1977 }
1978 }
1979 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr);
1980 }
1981 ierr = PetscSectionSetChart(parentSection,0,pNewEnd[dim]);CHKERRQ(ierr);
1982 offset = 0;
1983 numNewCones = 0;
1984 for (d = 0; d <= dim; d++) {
1985 PetscInt kStart, kEnd, k;
1986 PetscInt p;
1987 PetscInt size;
1988
1989 for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1990 /* skip cell 0 */
1991 if (p == cell) continue;
1992 /* old cones to new cones */
1993 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr);
1994 newConeSizes[offset++] = size;
1995 numNewCones += size;
1996 }
1997
1998 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr);
1999 for (k = kStart; k < kEnd; k++) {
2000 PetscInt kParent;
2001
2002 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr);
2003 if (kParent != k) {
2004 Kembedding[k] = offset;
2005 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr);
2006 newConeSizes[offset++] = size;
2007 numNewCones += size;
2008 if (kParent != 0) {
2009 ierr = PetscSectionSetDof(parentSection,Kembedding[k],1);CHKERRQ(ierr);
2010 }
2011 }
2012 }
2013 }
2014
2015 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr);
2016 ierr = PetscSectionGetStorageSize(parentSection,&numPointsWithParents);CHKERRQ(ierr);
2017 ierr = PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);CHKERRQ(ierr);
2018 ierr = PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);CHKERRQ(ierr);
2019
2020 /* fill new cones */
2021 offset = 0;
2022 for (d = 0; d <= dim; d++) {
2023 PetscInt kStart, kEnd, k, l;
2024 PetscInt p;
2025 PetscInt size;
2026 const PetscInt *cone, *orientation;
2027
2028 for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2029 /* skip cell 0 */
2030 if (p == cell) continue;
2031 /* old cones to new cones */
2032 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr);
2033 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr);
2034 ierr = DMPlexGetConeOrientation(dm,p,&orientation);CHKERRQ(ierr);
2035 for (l = 0; l < size; l++) {
2036 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2037 newOrientations[offset++] = orientation[l];
2038 }
2039 }
2040
2041 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr);
2042 for (k = kStart; k < kEnd; k++) {
2043 PetscInt kPerm = perm[k], kParent;
2044 PetscInt preO = preOrient[k];
2045
2046 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr);
2047 if (kParent != k) {
2048 /* embed new cones */
2049 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr);
2050 ierr = DMPlexGetCone(K,kPerm,&cone);CHKERRQ(ierr);
2051 ierr = DMPlexGetConeOrientation(K,kPerm,&orientation);CHKERRQ(ierr);
2052 for (l = 0; l < size; l++) {
2053 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2054 PetscInt newO, lSize, oTrue;
2055
2056 q = iperm[cone[m]];
2057 newCones[offset] = Kembedding[q];
2058 ierr = DMPlexGetConeSize(K,q,&lSize);CHKERRQ(ierr);
2059 oTrue = orientation[m];
2060 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2061 newO = DihedralCompose(lSize,oTrue,preOrient[q]);
2062 newOrientations[offset++] = newO;
2063 }
2064 if (kParent != 0) {
2065 PetscInt newPoint = Kembedding[kParent];
2066 ierr = PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);CHKERRQ(ierr);
2067 parents[pOffset] = newPoint;
2068 childIDs[pOffset] = k;
2069 }
2070 }
2071 }
2072 }
2073
2074 ierr = PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);CHKERRQ(ierr);
2075
2076 /* fill coordinates */
2077 offset = 0;
2078 {
2079 PetscInt kStart, kEnd, l;
2080 PetscSection vSection;
2081 PetscInt v;
2082 Vec coords;
2083 PetscScalar *coordvals;
2084 PetscInt dof, off;
2085 PetscReal v0[3], J[9], detJ;
2086
2087 if (PetscDefined(USE_DEBUG)) {
2088 PetscInt k;
2089 ierr = DMPlexGetHeightStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr);
2090 for (k = kStart; k < kEnd; k++) {
2091 ierr = DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr);
2092 if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2093 }
2094 }
2095 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr);
2096 ierr = DMGetCoordinateSection(dm,&vSection);CHKERRQ(ierr);
2097 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr);
2098 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr);
2099 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
2100
2101 ierr = PetscSectionGetDof(vSection,v,&dof);CHKERRQ(ierr);
2102 ierr = PetscSectionGetOffset(vSection,v,&off);CHKERRQ(ierr);
2103 for (l = 0; l < dof; l++) {
2104 newVertexCoords[offset++] = coordvals[off + l];
2105 }
2106 }
2107 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr);
2108
2109 ierr = DMGetCoordinateSection(K,&vSection);CHKERRQ(ierr);
2110 ierr = DMGetCoordinatesLocal(K,&coords);CHKERRQ(ierr);
2111 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr);
2112 ierr = DMPlexGetDepthStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr);
2113 for (v = kStart; v < kEnd; v++) {
2114 PetscReal coord[3], newCoord[3];
2115 PetscInt vPerm = perm[v];
2116 PetscInt kParent;
2117 const PetscReal xi0[3] = {-1.,-1.,-1.};
2118
2119 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr);
2120 if (kParent != v) {
2121 /* this is a new vertex */
2122 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr);
2123 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2124 CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2125 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2126 offset += dim;
2127 }
2128 }
2129 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr);
2130 }
2131
2132 /* need to reverse the order of pNewCount: vertices first, cells last */
2133 for (d = 0; d < (dim + 1) / 2; d++) {
2134 PetscInt tmp;
2135
2136 tmp = pNewCount[d];
2137 pNewCount[d] = pNewCount[dim - d];
2138 pNewCount[dim - d] = tmp;
2139 }
2140
2141 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr);
2142 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr);
2143 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr);
2144
2145 /* clean up */
2146 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr);
2147 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr);
2148 ierr = PetscFree(newConeSizes);CHKERRQ(ierr);
2149 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr);
2150 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr);
2151 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr);
2152 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr);
2153 }
2154 else {
2155 PetscInt p, counts[4];
2156 PetscInt *coneSizes, *cones, *orientations;
2157 Vec coordVec;
2158 PetscScalar *coords;
2159
2160 for (d = 0; d <= dim; d++) {
2161 PetscInt dStart, dEnd;
2162
2163 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr);
2164 counts[d] = dEnd - dStart;
2165 }
2166 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr);
2167 for (p = pStart; p < pEnd; p++) {
2168 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr);
2169 }
2170 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr);
2171 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr);
2172 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr);
2173 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr);
2174
2175 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr);
2176 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr);
2177 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr);
2178 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr);
2179 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr);
2180 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr);
2181 }
2182 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr);
2183
2184 PetscFunctionReturn(0);
2185 }
2186
DMPlexComputeInterpolatorTree(DM coarse,DM fine,PetscSF coarseToFine,PetscInt * childIds,Mat mat)2187 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2188 {
2189 PetscSF coarseToFineEmbedded;
2190 PetscSection globalCoarse, globalFine;
2191 PetscSection localCoarse, localFine;
2192 PetscSection aSec, cSec;
2193 PetscSection rootIndicesSec, rootMatricesSec;
2194 PetscSection leafIndicesSec, leafMatricesSec;
2195 PetscInt *rootIndices, *leafIndices;
2196 PetscScalar *rootMatrices, *leafMatrices;
2197 IS aIS;
2198 const PetscInt *anchors;
2199 Mat cMat;
2200 PetscInt numFields, maxFields;
2201 PetscInt pStartC, pEndC, pStartF, pEndF, p;
2202 PetscInt aStart, aEnd, cStart, cEnd;
2203 PetscInt *maxChildIds;
2204 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2205 const PetscInt ***perms;
2206 const PetscScalar ***flips;
2207 PetscErrorCode ierr;
2208
2209 PetscFunctionBegin;
2210 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
2211 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr);
2212 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr);
2213 { /* winnow fine points that don't have global dofs out of the sf */
2214 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2215 const PetscInt *leaves;
2216
2217 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr);
2218 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2219 p = leaves ? leaves[l] : l;
2220 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
2221 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
2222 if ((dof - cdof) > 0) {
2223 numPointsWithDofs++;
2224 }
2225 }
2226 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr);
2227 for (l = 0, offset = 0; l < nleaves; l++) {
2228 p = leaves ? leaves[l] : l;
2229 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
2230 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
2231 if ((dof - cdof) > 0) {
2232 pointsWithDofs[offset++] = l;
2233 }
2234 }
2235 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr);
2236 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr);
2237 }
2238 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2239 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr);
2240 for (p = pStartC; p < pEndC; p++) {
2241 maxChildIds[p - pStartC] = -2;
2242 }
2243 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr);
2244 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr);
2245
2246 ierr = DMGetLocalSection(coarse,&localCoarse);CHKERRQ(ierr);
2247 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr);
2248
2249 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr);
2250 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr);
2251 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
2252
2253 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr);
2254 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr);
2255
2256 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2257 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr);
2258 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr);
2259 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr);
2260 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr);
2261 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr);
2262 maxFields = PetscMax(1,numFields);
2263 ierr = PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);CHKERRQ(ierr);
2264 ierr = PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);CHKERRQ(ierr);
2265 ierr = PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));CHKERRQ(ierr);
2266 ierr = PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));CHKERRQ(ierr);
2267
2268 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2269 PetscInt dof, matSize = 0;
2270 PetscInt aDof = 0;
2271 PetscInt cDof = 0;
2272 PetscInt maxChildId = maxChildIds[p - pStartC];
2273 PetscInt numRowIndices = 0;
2274 PetscInt numColIndices = 0;
2275 PetscInt f;
2276
2277 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
2278 if (dof < 0) {
2279 dof = -(dof + 1);
2280 }
2281 if (p >= aStart && p < aEnd) {
2282 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr);
2283 }
2284 if (p >= cStart && p < cEnd) {
2285 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr);
2286 }
2287 for (f = 0; f <= numFields; f++) offsets[f] = 0;
2288 for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2289 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2290 PetscInt *closure = NULL, closureSize, cl;
2291
2292 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
2293 for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2294 PetscInt c = closure[2 * cl], clDof;
2295
2296 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr);
2297 numRowIndices += clDof;
2298 for (f = 0; f < numFields; f++) {
2299 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr);
2300 offsets[f + 1] += clDof;
2301 }
2302 }
2303 for (f = 0; f < numFields; f++) {
2304 offsets[f + 1] += offsets[f];
2305 newOffsets[f + 1] = offsets[f + 1];
2306 }
2307 /* get the number of indices needed and their field offsets */
2308 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr);
2309 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
2310 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2311 numColIndices = numRowIndices;
2312 matSize = 0;
2313 }
2314 else if (numFields) { /* we send one submat for each field: sum their sizes */
2315 matSize = 0;
2316 for (f = 0; f < numFields; f++) {
2317 PetscInt numRow, numCol;
2318
2319 numRow = offsets[f + 1] - offsets[f];
2320 numCol = newOffsets[f + 1] - newOffsets[f];
2321 matSize += numRow * numCol;
2322 }
2323 }
2324 else {
2325 matSize = numRowIndices * numColIndices;
2326 }
2327 } else if (maxChildId == -1) {
2328 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2329 PetscInt aOff, a;
2330
2331 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr);
2332 for (f = 0; f < numFields; f++) {
2333 PetscInt fDof;
2334
2335 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr);
2336 offsets[f+1] = fDof;
2337 }
2338 for (a = 0; a < aDof; a++) {
2339 PetscInt anchor = anchors[a + aOff], aLocalDof;
2340
2341 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr);
2342 numColIndices += aLocalDof;
2343 for (f = 0; f < numFields; f++) {
2344 PetscInt fDof;
2345
2346 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr);
2347 newOffsets[f+1] += fDof;
2348 }
2349 }
2350 if (numFields) {
2351 matSize = 0;
2352 for (f = 0; f < numFields; f++) {
2353 matSize += offsets[f+1] * newOffsets[f+1];
2354 }
2355 }
2356 else {
2357 matSize = numColIndices * dof;
2358 }
2359 }
2360 else { /* no children, and no constraints on dofs: just get the global indices */
2361 numColIndices = dof;
2362 matSize = 0;
2363 }
2364 }
2365 /* we will pack the column indices with the field offsets */
2366 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr);
2367 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr);
2368 }
2369 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr);
2370 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr);
2371 {
2372 PetscInt numRootIndices, numRootMatrices;
2373
2374 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr);
2375 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr);
2376 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr);
2377 for (p = pStartC; p < pEndC; p++) {
2378 PetscInt numRowIndices, numColIndices, matSize, dof;
2379 PetscInt pIndOff, pMatOff, f;
2380 PetscInt *pInd;
2381 PetscInt maxChildId = maxChildIds[p - pStartC];
2382 PetscScalar *pMat = NULL;
2383
2384 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr);
2385 if (!numColIndices) {
2386 continue;
2387 }
2388 for (f = 0; f <= numFields; f++) {
2389 offsets[f] = 0;
2390 newOffsets[f] = 0;
2391 offsetsCopy[f] = 0;
2392 newOffsetsCopy[f] = 0;
2393 }
2394 numColIndices -= 2 * numFields;
2395 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr);
2396 pInd = &(rootIndices[pIndOff]);
2397 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr);
2398 if (matSize) {
2399 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr);
2400 pMat = &rootMatrices[pMatOff];
2401 }
2402 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
2403 if (dof < 0) {
2404 dof = -(dof + 1);
2405 }
2406 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2407 PetscInt i, j;
2408 PetscInt numRowIndices = matSize / numColIndices;
2409
2410 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2411 PetscInt numIndices, *indices;
2412 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL);CHKERRQ(ierr);
2413 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2414 for (i = 0; i < numColIndices; i++) {
2415 pInd[i] = indices[i];
2416 }
2417 for (i = 0; i < numFields; i++) {
2418 pInd[numColIndices + i] = offsets[i+1];
2419 pInd[numColIndices + numFields + i] = offsets[i+1];
2420 }
2421 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL);CHKERRQ(ierr);
2422 }
2423 else {
2424 PetscInt closureSize, *closure = NULL, cl;
2425 PetscScalar *pMatIn, *pMatModified;
2426 PetscInt numPoints,*points;
2427
2428 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);CHKERRQ(ierr);
2429 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2430 for (j = 0; j < numRowIndices; j++) {
2431 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2432 }
2433 }
2434 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
2435 for (f = 0; f < maxFields; f++) {
2436 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);}
2437 else {ierr = PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);}
2438 }
2439 if (numFields) {
2440 for (cl = 0; cl < closureSize; cl++) {
2441 PetscInt c = closure[2 * cl];
2442
2443 for (f = 0; f < numFields; f++) {
2444 PetscInt fDof;
2445
2446 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr);
2447 offsets[f + 1] += fDof;
2448 }
2449 }
2450 for (f = 0; f < numFields; f++) {
2451 offsets[f + 1] += offsets[f];
2452 newOffsets[f + 1] = offsets[f + 1];
2453 }
2454 }
2455 /* TODO : flips here ? */
2456 /* apply hanging node constraints on the right, get the new points and the new offsets */
2457 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr);
2458 for (f = 0; f < maxFields; f++) {
2459 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);}
2460 else {ierr = PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);}
2461 }
2462 for (f = 0; f < maxFields; f++) {
2463 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);}
2464 else {ierr = PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);}
2465 }
2466 if (!numFields) {
2467 for (i = 0; i < numRowIndices * numColIndices; i++) {
2468 pMat[i] = pMatModified[i];
2469 }
2470 }
2471 else {
2472 PetscInt i, j, count;
2473 for (f = 0, count = 0; f < numFields; f++) {
2474 for (i = offsets[f]; i < offsets[f+1]; i++) {
2475 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2476 pMat[count] = pMatModified[i * numColIndices + j];
2477 }
2478 }
2479 }
2480 }
2481 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified);CHKERRQ(ierr);
2482 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
2483 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn);CHKERRQ(ierr);
2484 if (numFields) {
2485 for (f = 0; f < numFields; f++) {
2486 pInd[numColIndices + f] = offsets[f+1];
2487 pInd[numColIndices + numFields + f] = newOffsets[f+1];
2488 }
2489 for (cl = 0; cl < numPoints; cl++) {
2490 PetscInt globalOff, c = points[2*cl];
2491 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr);
2492 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd);CHKERRQ(ierr);
2493 }
2494 } else {
2495 for (cl = 0; cl < numPoints; cl++) {
2496 PetscInt c = points[2*cl], globalOff;
2497 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2498
2499 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr);
2500 ierr = DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd);CHKERRQ(ierr);
2501 }
2502 }
2503 for (f = 0; f < maxFields; f++) {
2504 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);}
2505 else {ierr = PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);}
2506 }
2507 ierr = DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);CHKERRQ(ierr);
2508 }
2509 }
2510 else if (matSize) {
2511 PetscInt cOff;
2512 PetscInt *rowIndices, *colIndices, a, aDof, aOff;
2513
2514 numRowIndices = matSize / numColIndices;
2515 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs");
2516 ierr = DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);CHKERRQ(ierr);
2517 ierr = DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);CHKERRQ(ierr);
2518 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr);
2519 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr);
2520 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr);
2521 if (numFields) {
2522 for (f = 0; f < numFields; f++) {
2523 PetscInt fDof;
2524
2525 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr);
2526 offsets[f + 1] = fDof;
2527 for (a = 0; a < aDof; a++) {
2528 PetscInt anchor = anchors[a + aOff];
2529 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr);
2530 newOffsets[f + 1] += fDof;
2531 }
2532 }
2533 for (f = 0; f < numFields; f++) {
2534 offsets[f + 1] += offsets[f];
2535 offsetsCopy[f + 1] = offsets[f + 1];
2536 newOffsets[f + 1] += newOffsets[f];
2537 newOffsetsCopy[f + 1] = newOffsets[f + 1];
2538 }
2539 ierr = DMPlexGetIndicesPointFields_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1, NULL,rowIndices);CHKERRQ(ierr);
2540 for (a = 0; a < aDof; a++) {
2541 PetscInt anchor = anchors[a + aOff], lOff;
2542 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr);
2543 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1, NULL,colIndices);CHKERRQ(ierr);
2544 }
2545 }
2546 else {
2547 ierr = DMPlexGetIndicesPoint_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL, NULL,rowIndices);CHKERRQ(ierr);
2548 for (a = 0; a < aDof; a++) {
2549 PetscInt anchor = anchors[a + aOff], lOff;
2550 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr);
2551 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL, NULL,colIndices);CHKERRQ(ierr);
2552 }
2553 }
2554 if (numFields) {
2555 PetscInt count, a;
2556
2557 for (f = 0, count = 0; f < numFields; f++) {
2558 PetscInt iSize = offsets[f + 1] - offsets[f];
2559 PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2560 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr);
2561 count += iSize * jSize;
2562 pInd[numColIndices + f] = offsets[f+1];
2563 pInd[numColIndices + numFields + f] = newOffsets[f+1];
2564 }
2565 for (a = 0; a < aDof; a++) {
2566 PetscInt anchor = anchors[a + aOff];
2567 PetscInt gOff;
2568 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr);
2569 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1, NULL,pInd);CHKERRQ(ierr);
2570 }
2571 }
2572 else {
2573 PetscInt a;
2574 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr);
2575 for (a = 0; a < aDof; a++) {
2576 PetscInt anchor = anchors[a + aOff];
2577 PetscInt gOff;
2578 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr);
2579 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL, NULL,pInd);CHKERRQ(ierr);
2580 }
2581 }
2582 ierr = DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);CHKERRQ(ierr);
2583 ierr = DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);CHKERRQ(ierr);
2584 }
2585 else {
2586 PetscInt gOff;
2587
2588 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr);
2589 if (numFields) {
2590 for (f = 0; f < numFields; f++) {
2591 PetscInt fDof;
2592 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr);
2593 offsets[f + 1] = fDof + offsets[f];
2594 }
2595 for (f = 0; f < numFields; f++) {
2596 pInd[numColIndices + f] = offsets[f+1];
2597 pInd[numColIndices + numFields + f] = offsets[f+1];
2598 }
2599 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd);CHKERRQ(ierr);
2600 } else {
2601 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd);CHKERRQ(ierr);
2602 }
2603 }
2604 }
2605 ierr = PetscFree(maxChildIds);CHKERRQ(ierr);
2606 }
2607 {
2608 PetscSF indicesSF, matricesSF;
2609 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2610
2611 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr);
2612 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr);
2613 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr);
2614 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr);
2615 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr);
2616 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr);
2617 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr);
2618 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr);
2619 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr);
2620 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr);
2621 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr);
2622 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr);
2623 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr);
2624 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr);
2625 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr);
2626 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr);
2627 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr);
2628 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr);
2629 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr);
2630 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr);
2631 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr);
2632 }
2633 /* count to preallocate */
2634 ierr = DMGetLocalSection(fine,&localFine);CHKERRQ(ierr);
2635 {
2636 PetscInt nGlobal;
2637 PetscInt *dnnz, *onnz;
2638 PetscLayout rowMap, colMap;
2639 PetscInt rowStart, rowEnd, colStart, colEnd;
2640 PetscInt maxDof;
2641 PetscInt *rowIndices;
2642 DM refTree;
2643 PetscInt **refPointFieldN;
2644 PetscScalar ***refPointFieldMats;
2645 PetscSection refConSec, refAnSec;
2646 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd;
2647 PetscScalar *pointWork;
2648
2649 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr);
2650 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr);
2651 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr);
2652 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr);
2653 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr);
2654 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr);
2655 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr);
2656 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr);
2657 ierr = PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);CHKERRQ(ierr);
2658 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr);
2659 for (p = leafStart; p < leafEnd; p++) {
2660 PetscInt gDof, gcDof, gOff;
2661 PetscInt numColIndices, pIndOff, *pInd;
2662 PetscInt matSize;
2663 PetscInt i;
2664
2665 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr);
2666 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr);
2667 if ((gDof - gcDof) <= 0) {
2668 continue;
2669 }
2670 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr);
2671 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset");
2672 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs");
2673 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr);
2674 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr);
2675 numColIndices -= 2 * numFields;
2676 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from");
2677 pInd = &leafIndices[pIndOff];
2678 offsets[0] = 0;
2679 offsetsCopy[0] = 0;
2680 newOffsets[0] = 0;
2681 newOffsetsCopy[0] = 0;
2682 if (numFields) {
2683 PetscInt f;
2684 for (f = 0; f < numFields; f++) {
2685 PetscInt rowDof;
2686
2687 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr);
2688 offsets[f + 1] = offsets[f] + rowDof;
2689 offsetsCopy[f + 1] = offsets[f + 1];
2690 newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2691 numD[f] = 0;
2692 numO[f] = 0;
2693 }
2694 ierr = DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices);CHKERRQ(ierr);
2695 for (f = 0; f < numFields; f++) {
2696 PetscInt colOffset = newOffsets[f];
2697 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2698
2699 for (i = 0; i < numFieldCols; i++) {
2700 PetscInt gInd = pInd[i + colOffset];
2701
2702 if (gInd >= colStart && gInd < colEnd) {
2703 numD[f]++;
2704 }
2705 else if (gInd >= 0) { /* negative means non-entry */
2706 numO[f]++;
2707 }
2708 }
2709 }
2710 }
2711 else {
2712 ierr = DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices);CHKERRQ(ierr);
2713 numD[0] = 0;
2714 numO[0] = 0;
2715 for (i = 0; i < numColIndices; i++) {
2716 PetscInt gInd = pInd[i];
2717
2718 if (gInd >= colStart && gInd < colEnd) {
2719 numD[0]++;
2720 }
2721 else if (gInd >= 0) { /* negative means non-entry */
2722 numO[0]++;
2723 }
2724 }
2725 }
2726 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr);
2727 if (!matSize) { /* incoming matrix is identity */
2728 PetscInt childId;
2729
2730 childId = childIds[p-pStartF];
2731 if (childId < 0) { /* no child interpolation: one nnz per */
2732 if (numFields) {
2733 PetscInt f;
2734 for (f = 0; f < numFields; f++) {
2735 PetscInt numRows = offsets[f+1] - offsets[f], row;
2736 for (row = 0; row < numRows; row++) {
2737 PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2738 PetscInt gIndFine = rowIndices[offsets[f] + row];
2739 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2740 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2741 dnnz[gIndFine - rowStart] = 1;
2742 }
2743 else if (gIndCoarse >= 0) { /* remote */
2744 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2745 onnz[gIndFine - rowStart] = 1;
2746 }
2747 else { /* constrained */
2748 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2749 }
2750 }
2751 }
2752 }
2753 else {
2754 PetscInt i;
2755 for (i = 0; i < gDof; i++) {
2756 PetscInt gIndCoarse = pInd[i];
2757 PetscInt gIndFine = rowIndices[i];
2758 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2759 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2760 dnnz[gIndFine - rowStart] = 1;
2761 }
2762 else if (gIndCoarse >= 0) { /* remote */
2763 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2764 onnz[gIndFine - rowStart] = 1;
2765 }
2766 else { /* constrained */
2767 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2768 }
2769 }
2770 }
2771 }
2772 else { /* interpolate from all */
2773 if (numFields) {
2774 PetscInt f;
2775 for (f = 0; f < numFields; f++) {
2776 PetscInt numRows = offsets[f+1] - offsets[f], row;
2777 for (row = 0; row < numRows; row++) {
2778 PetscInt gIndFine = rowIndices[offsets[f] + row];
2779 if (gIndFine >= 0) {
2780 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2781 dnnz[gIndFine - rowStart] = numD[f];
2782 onnz[gIndFine - rowStart] = numO[f];
2783 }
2784 }
2785 }
2786 }
2787 else {
2788 PetscInt i;
2789 for (i = 0; i < gDof; i++) {
2790 PetscInt gIndFine = rowIndices[i];
2791 if (gIndFine >= 0) {
2792 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2793 dnnz[gIndFine - rowStart] = numD[0];
2794 onnz[gIndFine - rowStart] = numO[0];
2795 }
2796 }
2797 }
2798 }
2799 }
2800 else { /* interpolate from all */
2801 if (numFields) {
2802 PetscInt f;
2803 for (f = 0; f < numFields; f++) {
2804 PetscInt numRows = offsets[f+1] - offsets[f], row;
2805 for (row = 0; row < numRows; row++) {
2806 PetscInt gIndFine = rowIndices[offsets[f] + row];
2807 if (gIndFine >= 0) {
2808 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2809 dnnz[gIndFine - rowStart] = numD[f];
2810 onnz[gIndFine - rowStart] = numO[f];
2811 }
2812 }
2813 }
2814 }
2815 else { /* every dof get a full row */
2816 PetscInt i;
2817 for (i = 0; i < gDof; i++) {
2818 PetscInt gIndFine = rowIndices[i];
2819 if (gIndFine >= 0) {
2820 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2821 dnnz[gIndFine - rowStart] = numD[0];
2822 onnz[gIndFine - rowStart] = numO[0];
2823 }
2824 }
2825 }
2826 }
2827 }
2828 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr);
2829 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr);
2830
2831 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr);
2832 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
2833 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr);
2834 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr);
2835 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
2836 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr);
2837 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr);
2838 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr);
2839 for (p = leafStart; p < leafEnd; p++) {
2840 PetscInt gDof, gcDof, gOff;
2841 PetscInt numColIndices, pIndOff, *pInd;
2842 PetscInt matSize;
2843 PetscInt childId;
2844
2845 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr);
2846 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr);
2847 if ((gDof - gcDof) <= 0) {
2848 continue;
2849 }
2850 childId = childIds[p-pStartF];
2851 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr);
2852 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr);
2853 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr);
2854 numColIndices -= 2 * numFields;
2855 pInd = &leafIndices[pIndOff];
2856 offsets[0] = 0;
2857 offsetsCopy[0] = 0;
2858 newOffsets[0] = 0;
2859 newOffsetsCopy[0] = 0;
2860 rowOffsets[0] = 0;
2861 if (numFields) {
2862 PetscInt f;
2863 for (f = 0; f < numFields; f++) {
2864 PetscInt rowDof;
2865
2866 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr);
2867 offsets[f + 1] = offsets[f] + rowDof;
2868 offsetsCopy[f + 1] = offsets[f + 1];
2869 rowOffsets[f + 1] = pInd[numColIndices + f];
2870 newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2871 }
2872 ierr = DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices);CHKERRQ(ierr);
2873 }
2874 else {
2875 ierr = DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices);CHKERRQ(ierr);
2876 }
2877 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr);
2878 if (!matSize) { /* incoming matrix is identity */
2879 if (childId < 0) { /* no child interpolation: scatter */
2880 if (numFields) {
2881 PetscInt f;
2882 for (f = 0; f < numFields; f++) {
2883 PetscInt numRows = offsets[f+1] - offsets[f], row;
2884 for (row = 0; row < numRows; row++) {
2885 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr);
2886 }
2887 }
2888 }
2889 else {
2890 PetscInt numRows = gDof, row;
2891 for (row = 0; row < numRows; row++) {
2892 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr);
2893 }
2894 }
2895 }
2896 else { /* interpolate from all */
2897 if (numFields) {
2898 PetscInt f;
2899 for (f = 0; f < numFields; f++) {
2900 PetscInt numRows = offsets[f+1] - offsets[f];
2901 PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2902 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr);
2903 }
2904 }
2905 else {
2906 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr);
2907 }
2908 }
2909 }
2910 else { /* interpolate from all */
2911 PetscInt pMatOff;
2912 PetscScalar *pMat;
2913
2914 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr);
2915 pMat = &leafMatrices[pMatOff];
2916 if (childId < 0) { /* copy the incoming matrix */
2917 if (numFields) {
2918 PetscInt f, count;
2919 for (f = 0, count = 0; f < numFields; f++) {
2920 PetscInt numRows = offsets[f+1]-offsets[f];
2921 PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2922 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2923 PetscScalar *inMat = &pMat[count];
2924
2925 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr);
2926 count += numCols * numInRows;
2927 }
2928 }
2929 else {
2930 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr);
2931 }
2932 }
2933 else { /* multiply the incoming matrix by the child interpolation */
2934 if (numFields) {
2935 PetscInt f, count;
2936 for (f = 0, count = 0; f < numFields; f++) {
2937 PetscInt numRows = offsets[f+1]-offsets[f];
2938 PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2939 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2940 PetscScalar *inMat = &pMat[count];
2941 PetscInt i, j, k;
2942 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2943 for (i = 0; i < numRows; i++) {
2944 for (j = 0; j < numCols; j++) {
2945 PetscScalar val = 0.;
2946 for (k = 0; k < numInRows; k++) {
2947 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2948 }
2949 pointWork[i * numCols + j] = val;
2950 }
2951 }
2952 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr);
2953 count += numCols * numInRows;
2954 }
2955 }
2956 else { /* every dof gets a full row */
2957 PetscInt numRows = gDof;
2958 PetscInt numCols = numColIndices;
2959 PetscInt numInRows = matSize / numColIndices;
2960 PetscInt i, j, k;
2961 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2962 for (i = 0; i < numRows; i++) {
2963 for (j = 0; j < numCols; j++) {
2964 PetscScalar val = 0.;
2965 for (k = 0; k < numInRows; k++) {
2966 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2967 }
2968 pointWork[i * numCols + j] = val;
2969 }
2970 }
2971 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr);
2972 }
2973 }
2974 }
2975 }
2976 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
2977 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr);
2978 ierr = PetscFree(pointWork);CHKERRQ(ierr);
2979 }
2980 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2981 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2982 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr);
2983 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr);
2984 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr);
2985 ierr = PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);CHKERRQ(ierr);
2986 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr);
2987 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr);
2988 PetscFunctionReturn(0);
2989 }
2990
2991 /*
2992 * Assuming a nodal basis (w.r.t. the dual basis) basis:
2993 *
2994 * for each coarse dof \phi^c_i:
2995 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2996 * for each fine dof \phi^f_j;
2997 * a_{i,j} = 0;
2998 * for each fine dof \phi^f_k:
2999 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
3000 * [^^^ this is = \phi^c_i ^^^]
3001 */
DMPlexComputeInjectorReferenceTree(DM refTree,Mat * inj)3002 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
3003 {
3004 PetscDS ds;
3005 PetscSection section, cSection;
3006 DMLabel canonical, depth;
3007 Mat cMat, mat;
3008 PetscInt *nnz;
3009 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
3010 PetscInt m, n;
3011 PetscScalar *pointScalar;
3012 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
3013 PetscErrorCode ierr;
3014
3015 PetscFunctionBegin;
3016 ierr = DMGetLocalSection(refTree,§ion);CHKERRQ(ierr);
3017 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr);
3018 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr);
3019 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr);
3020 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr);
3021 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
3022 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr);
3023 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr);
3024 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr);
3025 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr);
3026 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr);
3027 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr);
3028 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */
3029 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
3030 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr);
3031 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
3032 const PetscInt *children;
3033 PetscInt numChildren;
3034 PetscInt i, numChildDof, numSelfDof;
3035
3036 if (canonical) {
3037 PetscInt pCanonical;
3038 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr);
3039 if (p != pCanonical) continue;
3040 }
3041 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr);
3042 if (!numChildren) continue;
3043 for (i = 0, numChildDof = 0; i < numChildren; i++) {
3044 PetscInt child = children[i];
3045 PetscInt dof;
3046
3047 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr);
3048 numChildDof += dof;
3049 }
3050 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr);
3051 if (!numChildDof || !numSelfDof) continue;
3052 for (f = 0; f < numFields; f++) {
3053 PetscInt selfOff;
3054
3055 if (numSecFields) { /* count the dofs for just this field */
3056 for (i = 0, numChildDof = 0; i < numChildren; i++) {
3057 PetscInt child = children[i];
3058 PetscInt dof;
3059
3060 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr);
3061 numChildDof += dof;
3062 }
3063 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr);
3064 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr);
3065 }
3066 else {
3067 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr);
3068 }
3069 for (i = 0; i < numSelfDof; i++) {
3070 nnz[selfOff + i] = numChildDof;
3071 }
3072 }
3073 }
3074 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr);
3075 ierr = PetscFree(nnz);CHKERRQ(ierr);
3076 /* Setp 2: compute entries */
3077 for (p = pStart; p < pEnd; p++) {
3078 const PetscInt *children;
3079 PetscInt numChildren;
3080 PetscInt i, numChildDof, numSelfDof;
3081
3082 /* same conditions about when entries occur */
3083 if (canonical) {
3084 PetscInt pCanonical;
3085 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr);
3086 if (p != pCanonical) continue;
3087 }
3088 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr);
3089 if (!numChildren) continue;
3090 for (i = 0, numChildDof = 0; i < numChildren; i++) {
3091 PetscInt child = children[i];
3092 PetscInt dof;
3093
3094 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr);
3095 numChildDof += dof;
3096 }
3097 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr);
3098 if (!numChildDof || !numSelfDof) continue;
3099
3100 for (f = 0; f < numFields; f++) {
3101 PetscInt pI = -1, cI = -1;
3102 PetscInt selfOff, Nc, parentCell;
3103 PetscInt cellShapeOff;
3104 PetscObject disc;
3105 PetscDualSpace dsp;
3106 PetscClassId classId;
3107 PetscScalar *pointMat;
3108 PetscInt *matRows, *matCols;
3109 PetscInt pO = PETSC_MIN_INT;
3110 const PetscInt *depthNumDof;
3111
3112 if (numSecFields) {
3113 for (i = 0, numChildDof = 0; i < numChildren; i++) {
3114 PetscInt child = children[i];
3115 PetscInt dof;
3116
3117 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr);
3118 numChildDof += dof;
3119 }
3120 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr);
3121 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr);
3122 }
3123 else {
3124 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr);
3125 }
3126
3127 /* find a cell whose closure contains p */
3128 if (p >= cStart && p < cEnd) {
3129 parentCell = p;
3130 }
3131 else {
3132 PetscInt *star = NULL;
3133 PetscInt numStar;
3134
3135 parentCell = -1;
3136 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr);
3137 for (i = numStar - 1; i >= 0; i--) {
3138 PetscInt c = star[2 * i];
3139
3140 if (c >= cStart && c < cEnd) {
3141 parentCell = c;
3142 break;
3143 }
3144 }
3145 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr);
3146 }
3147 /* determine the offset of p's shape functions withing parentCell's shape functions */
3148 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr);
3149 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr);
3150 if (classId == PETSCFE_CLASSID) {
3151 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr);
3152 }
3153 else if (classId == PETSCFV_CLASSID) {
3154 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr);
3155 }
3156 else {
3157 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3158 }
3159 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr);
3160 ierr = PetscDualSpaceGetNumComponents(dsp,&Nc);CHKERRQ(ierr);
3161 {
3162 PetscInt *closure = NULL;
3163 PetscInt numClosure;
3164
3165 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr);
3166 for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
3167 PetscInt point = closure[2 * i], pointDepth;
3168
3169 pO = closure[2 * i + 1];
3170 if (point == p) {
3171 pI = i;
3172 break;
3173 }
3174 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr);
3175 cellShapeOff += depthNumDof[pointDepth];
3176 }
3177 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr);
3178 }
3179
3180 ierr = DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr);
3181 ierr = DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr);
3182 matCols = matRows + numSelfDof;
3183 for (i = 0; i < numSelfDof; i++) {
3184 matRows[i] = selfOff + i;
3185 }
3186 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3187 {
3188 PetscInt colOff = 0;
3189
3190 for (i = 0; i < numChildren; i++) {
3191 PetscInt child = children[i];
3192 PetscInt dof, off, j;
3193
3194 if (numSecFields) {
3195 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr);
3196 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr);
3197 }
3198 else {
3199 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr);
3200 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr);
3201 }
3202
3203 for (j = 0; j < dof; j++) {
3204 matCols[colOff++] = off + j;
3205 }
3206 }
3207 }
3208 if (classId == PETSCFE_CLASSID) {
3209 PetscFE fe = (PetscFE) disc;
3210 PetscInt fSize;
3211 const PetscInt ***perms;
3212 const PetscScalar ***flips;
3213 const PetscInt *pperms;
3214
3215
3216 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr);
3217 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr);
3218 ierr = PetscDualSpaceGetSymmetries(dsp, &perms, &flips);CHKERRQ(ierr);
3219 pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3220 for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3221 PetscQuadrature q;
3222 PetscInt dim, thisNc, numPoints, j, k;
3223 const PetscReal *points;
3224 const PetscReal *weights;
3225 PetscInt *closure = NULL;
3226 PetscInt numClosure;
3227 PetscInt iCell = pperms ? pperms[i] : i;
3228 PetscInt parentCellShapeDof = cellShapeOff + iCell;
3229 PetscTabulation Tparent;
3230
3231 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr);
3232 ierr = PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);CHKERRQ(ierr);
3233 if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
3234 ierr = PetscFECreateTabulation(fe,1,numPoints,points,0,&Tparent);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3235 for (j = 0; j < numPoints; j++) {
3236 PetscInt childCell = -1;
3237 PetscReal *parentValAtPoint;
3238 const PetscReal xi0[3] = {-1.,-1.,-1.};
3239 const PetscReal *pointReal = &points[dim * j];
3240 const PetscScalar *point;
3241 PetscTabulation Tchild;
3242 PetscInt childCellShapeOff, pointMatOff;
3243 #if defined(PETSC_USE_COMPLEX)
3244 PetscInt d;
3245
3246 for (d = 0; d < dim; d++) {
3247 pointScalar[d] = points[dim * j + d];
3248 }
3249 point = pointScalar;
3250 #else
3251 point = pointReal;
3252 #endif
3253
3254 parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];
3255
3256 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3257 PetscInt child = children[k];
3258 PetscInt *star = NULL;
3259 PetscInt numStar, s;
3260
3261 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr);
3262 for (s = numStar - 1; s >= 0; s--) {
3263 PetscInt c = star[2 * s];
3264
3265 if (c < cStart || c >= cEnd) continue;
3266 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr);
3267 if (childCell >= 0) break;
3268 }
3269 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr);
3270 if (childCell >= 0) break;
3271 }
3272 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");
3273 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
3274 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr);
3275 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3276 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3277
3278 ierr = PetscFECreateTabulation(fe,1,1,pointRef,0,&Tchild);CHKERRQ(ierr);
3279 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr);
3280 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3281 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3282 PetscInt l;
3283 const PetscInt *cperms;
3284
3285 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr);
3286 childDof = depthNumDof[childDepth];
3287 for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3288 PetscInt point = closure[2 * l];
3289 PetscInt pointDepth;
3290
3291 childO = closure[2 * l + 1];
3292 if (point == child) {
3293 cI = l;
3294 break;
3295 }
3296 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr);
3297 childCellShapeOff += depthNumDof[pointDepth];
3298 }
3299 if (l == numClosure) {
3300 pointMatOff += childDof;
3301 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3302 }
3303 cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3304 for (l = 0; l < childDof; l++) {
3305 PetscInt lCell = cperms ? cperms[l] : l;
3306 PetscInt childCellDof = childCellShapeOff + lCell;
3307 PetscReal *childValAtPoint;
3308 PetscReal val = 0.;
3309
3310 childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3311 for (m = 0; m < Nc; m++) {
3312 val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3313 }
3314
3315 pointMat[i * numChildDof + pointMatOff + l] += val;
3316 }
3317 pointMatOff += childDof;
3318 }
3319 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr);
3320 ierr = PetscTabulationDestroy(&Tchild);CHKERRQ(ierr);
3321 }
3322 ierr = PetscTabulationDestroy(&Tparent);CHKERRQ(ierr);
3323 }
3324 }
3325 else { /* just the volume-weighted averages of the children */
3326 PetscReal parentVol;
3327 PetscInt childCell;
3328
3329 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr);
3330 for (i = 0, childCell = 0; i < numChildren; i++) {
3331 PetscInt child = children[i], j;
3332 PetscReal childVol;
3333
3334 if (child < cStart || child >= cEnd) continue;
3335 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr);
3336 for (j = 0; j < Nc; j++) {
3337 pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3338 }
3339 childCell++;
3340 }
3341 }
3342 /* Insert pointMat into mat */
3343 ierr = MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr);
3344 ierr = DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr);
3345 ierr = DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr);
3346 }
3347 }
3348 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr);
3349 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr);
3350 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3351 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3352 *inj = mat;
3353 PetscFunctionReturn(0);
3354 }
3355
DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree,Mat inj,PetscScalar **** childrenMats)3356 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3357 {
3358 PetscDS ds;
3359 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3360 PetscScalar ***refPointFieldMats;
3361 PetscSection refConSec, refSection;
3362 PetscErrorCode ierr;
3363
3364 PetscFunctionBegin;
3365 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr);
3366 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
3367 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr);
3368 ierr = DMGetLocalSection(refTree,&refSection);CHKERRQ(ierr);
3369 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
3370 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr);
3371 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr);
3372 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr);
3373 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr);
3374 for (p = pRefStart; p < pRefEnd; p++) {
3375 PetscInt parent, pDof, parentDof;
3376
3377 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr);
3378 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr);
3379 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr);
3380 if (!pDof || !parentDof || parent == p) continue;
3381
3382 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr);
3383 for (f = 0; f < numFields; f++) {
3384 PetscInt cDof, cOff, numCols, r;
3385
3386 if (numFields > 1) {
3387 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr);
3388 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr);
3389 }
3390 else {
3391 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr);
3392 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr);
3393 }
3394
3395 for (r = 0; r < cDof; r++) {
3396 rows[r] = cOff + r;
3397 }
3398 numCols = 0;
3399 {
3400 PetscInt aDof, aOff, j;
3401
3402 if (numFields > 1) {
3403 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr);
3404 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr);
3405 }
3406 else {
3407 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr);
3408 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr);
3409 }
3410
3411 for (j = 0; j < aDof; j++) {
3412 cols[numCols++] = aOff + j;
3413 }
3414 }
3415 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr);
3416 /* transpose of constraint matrix */
3417 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr);
3418 }
3419 }
3420 *childrenMats = refPointFieldMats;
3421 ierr = PetscFree(rows);CHKERRQ(ierr);
3422 ierr = PetscFree(cols);CHKERRQ(ierr);
3423 PetscFunctionReturn(0);
3424 }
3425
DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree,Mat inj,PetscScalar **** childrenMats)3426 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3427 {
3428 PetscDS ds;
3429 PetscScalar ***refPointFieldMats;
3430 PetscInt numFields, pRefStart, pRefEnd, p, f;
3431 PetscSection refConSec, refSection;
3432 PetscErrorCode ierr;
3433
3434 PetscFunctionBegin;
3435 refPointFieldMats = *childrenMats;
3436 *childrenMats = NULL;
3437 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr);
3438 ierr = DMGetLocalSection(refTree,&refSection);CHKERRQ(ierr);
3439 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr);
3440 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr);
3441 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
3442 for (p = pRefStart; p < pRefEnd; p++) {
3443 PetscInt parent, pDof, parentDof;
3444
3445 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr);
3446 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr);
3447 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr);
3448 if (!pDof || !parentDof || parent == p) continue;
3449
3450 for (f = 0; f < numFields; f++) {
3451 PetscInt cDof;
3452
3453 if (numFields > 1) {
3454 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr);
3455 }
3456 else {
3457 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr);
3458 }
3459
3460 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr);
3461 }
3462 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr);
3463 }
3464 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr);
3465 PetscFunctionReturn(0);
3466 }
3467
DMPlexReferenceTreeGetInjector(DM refTree,Mat * injRef)3468 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef)
3469 {
3470 Mat cMatRef;
3471 PetscObject injRefObj;
3472 PetscErrorCode ierr;
3473
3474 PetscFunctionBegin;
3475 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr);
3476 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr);
3477 *injRef = (Mat) injRefObj;
3478 if (!*injRef) {
3479 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr);
3480 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr);
3481 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3482 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr);
3483 }
3484 PetscFunctionReturn(0);
3485 }
3486
DMPlexTransferInjectorTree(DM coarse,DM fine,PetscSF coarseToFine,const PetscInt * childIds,Vec fineVec,PetscInt numFields,PetscInt * offsets,PetscSection * rootMultiSec,PetscSection * multiLeafSec,PetscInt ** gatheredIndices,PetscScalar ** gatheredValues)3487 static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues)
3488 {
3489 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3490 PetscSection globalCoarse, globalFine;
3491 PetscSection localCoarse, localFine, leafIndicesSec;
3492 PetscSection multiRootSec, rootIndicesSec;
3493 PetscInt *leafInds, *rootInds = NULL;
3494 const PetscInt *rootDegrees;
3495 PetscScalar *leafVals = NULL, *rootVals = NULL;
3496 PetscSF coarseToFineEmbedded;
3497 PetscErrorCode ierr;
3498
3499 PetscFunctionBegin;
3500 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
3501 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr);
3502 ierr = DMGetLocalSection(fine,&localFine);CHKERRQ(ierr);
3503 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr);
3504 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr);
3505 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr);
3506 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr);
3507 { /* winnow fine points that don't have global dofs out of the sf */
3508 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3509 const PetscInt *leaves;
3510
3511 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr);
3512 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3513 p = leaves ? leaves[l] : l;
3514 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
3515 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
3516 if ((dof - cdof) > 0) {
3517 numPointsWithDofs++;
3518
3519 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr);
3520 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr);
3521 }
3522 }
3523 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr);
3524 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr);
3525 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr);
3526 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr);
3527 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);}
3528 for (l = 0, offset = 0; l < nleaves; l++) {
3529 p = leaves ? leaves[l] : l;
3530 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
3531 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
3532 if ((dof - cdof) > 0) {
3533 PetscInt off, gOff;
3534 PetscInt *pInd;
3535 PetscScalar *pVal = NULL;
3536
3537 pointsWithDofs[offset++] = l;
3538
3539 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr);
3540
3541 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3542 if (gatheredValues) {
3543 PetscInt i;
3544
3545 pVal = &leafVals[off + 1];
3546 for (i = 0; i < dof; i++) pVal[i] = 0.;
3547 }
3548 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr);
3549
3550 offsets[0] = 0;
3551 if (numFields) {
3552 PetscInt f;
3553
3554 for (f = 0; f < numFields; f++) {
3555 PetscInt fDof;
3556 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr);
3557 offsets[f + 1] = fDof + offsets[f];
3558 }
3559 ierr = DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd);CHKERRQ(ierr);
3560 } else {
3561 ierr = DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd);CHKERRQ(ierr);
3562 }
3563 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);}
3564 }
3565 }
3566 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr);
3567 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr);
3568 }
3569
3570 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
3571 ierr = DMGetLocalSection(coarse,&localCoarse);CHKERRQ(ierr);
3572 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr);
3573
3574 { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3575 MPI_Datatype threeInt;
3576 PetscMPIInt rank;
3577 PetscInt (*parentNodeAndIdCoarse)[3];
3578 PetscInt (*parentNodeAndIdFine)[3];
3579 PetscInt p, nleaves, nleavesToParents;
3580 PetscSF pointSF, sfToParents;
3581 const PetscInt *ilocal;
3582 const PetscSFNode *iremote;
3583 PetscSFNode *iremoteToParents;
3584 PetscInt *ilocalToParents;
3585
3586 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr);
3587 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr);
3588 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr);
3589 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr);
3590 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr);
3591 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr);
3592 for (p = pStartC; p < pEndC; p++) {
3593 PetscInt parent, childId;
3594 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr);
3595 parentNodeAndIdCoarse[p - pStartC][0] = rank;
3596 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3597 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3598 if (nleaves > 0) {
3599 PetscInt leaf = -1;
3600
3601 if (ilocal) {
3602 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr);
3603 }
3604 else {
3605 leaf = p - pStartC;
3606 }
3607 if (leaf >= 0) {
3608 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3609 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3610 }
3611 }
3612 }
3613 for (p = pStartF; p < pEndF; p++) {
3614 parentNodeAndIdFine[p - pStartF][0] = -1;
3615 parentNodeAndIdFine[p - pStartF][1] = -1;
3616 parentNodeAndIdFine[p - pStartF][2] = -1;
3617 }
3618 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr);
3619 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr);
3620 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3621 PetscInt dof;
3622
3623 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr);
3624 if (dof) {
3625 PetscInt off;
3626
3627 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr);
3628 if (gatheredIndices) {
3629 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3630 } else if (gatheredValues) {
3631 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3632 }
3633 }
3634 if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3635 nleavesToParents++;
3636 }
3637 }
3638 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr);
3639 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr);
3640 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3641 if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3642 ilocalToParents[nleavesToParents] = p - pStartF;
3643 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0];
3644 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3645 nleavesToParents++;
3646 }
3647 }
3648 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr);
3649 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr);
3650 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr);
3651
3652 coarseToFineEmbedded = sfToParents;
3653
3654 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr);
3655 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr);
3656 }
3657
3658 { /* winnow out coarse points that don't have dofs */
3659 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3660 PetscSF sfDofsOnly;
3661
3662 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3663 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
3664 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr);
3665 if ((dof - cdof) > 0) {
3666 numPointsWithDofs++;
3667 }
3668 }
3669 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr);
3670 for (p = pStartC, offset = 0; p < pEndC; p++) {
3671 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
3672 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr);
3673 if ((dof - cdof) > 0) {
3674 pointsWithDofs[offset++] = p - pStartC;
3675 }
3676 }
3677 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr);
3678 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr);
3679 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr);
3680 coarseToFineEmbedded = sfDofsOnly;
3681 }
3682
3683 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3684 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr);
3685 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr);
3686 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr);
3687 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr);
3688 for (p = pStartC; p < pEndC; p++) {
3689 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr);
3690 }
3691 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr);
3692 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr);
3693 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr);
3694 { /* distribute the leaf section */
3695 PetscSF multi, multiInv, indicesSF;
3696 PetscInt *remoteOffsets, numRootIndices;
3697
3698 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr);
3699 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr);
3700 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr);
3701 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr);
3702 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
3703 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr);
3704 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr);
3705 if (gatheredIndices) {
3706 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr);
3707 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr);
3708 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr);
3709 }
3710 if (gatheredValues) {
3711 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr);
3712 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr);
3713 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr);
3714 }
3715 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr);
3716 }
3717 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr);
3718 ierr = PetscFree(leafInds);CHKERRQ(ierr);
3719 ierr = PetscFree(leafVals);CHKERRQ(ierr);
3720 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr);
3721 *rootMultiSec = multiRootSec;
3722 *multiLeafSec = rootIndicesSec;
3723 if (gatheredIndices) *gatheredIndices = rootInds;
3724 if (gatheredValues) *gatheredValues = rootVals;
3725 PetscFunctionReturn(0);
3726 }
3727
DMPlexComputeInjectorTree(DM coarse,DM fine,PetscSF coarseToFine,PetscInt * childIds,Mat mat)3728 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3729 {
3730 DM refTree;
3731 PetscSection multiRootSec, rootIndicesSec;
3732 PetscSection globalCoarse, globalFine;
3733 PetscSection localCoarse, localFine;
3734 PetscSection cSecRef;
3735 PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3736 Mat injRef;
3737 PetscInt numFields, maxDof;
3738 PetscInt pStartC, pEndC, pStartF, pEndF, p;
3739 PetscInt *offsets, *offsetsCopy, *rowOffsets;
3740 PetscLayout rowMap, colMap;
3741 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3742 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3743 PetscErrorCode ierr;
3744
3745 PetscFunctionBegin;
3746
3747 /* get the templates for the fine-to-coarse injection from the reference tree */
3748 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr);
3749 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr);
3750 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr);
3751 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr);
3752
3753 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr);
3754 ierr = DMGetLocalSection(fine,&localFine);CHKERRQ(ierr);
3755 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr);
3756 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr);
3757 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
3758 ierr = DMGetLocalSection(coarse,&localCoarse);CHKERRQ(ierr);
3759 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr);
3760 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr);
3761 {
3762 PetscInt maxFields = PetscMax(1,numFields) + 1;
3763 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr);
3764 }
3765
3766 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr);
3767
3768 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr);
3769
3770 /* count indices */
3771 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr);
3772 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr);
3773 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr);
3774 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr);
3775 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr);
3776 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr);
3777 for (p = pStartC; p < pEndC; p++) {
3778 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3779
3780 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
3781 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr);
3782 if ((dof - cdof) <= 0) continue;
3783 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr);
3784
3785 rowOffsets[0] = 0;
3786 offsetsCopy[0] = 0;
3787 if (numFields) {
3788 PetscInt f;
3789
3790 for (f = 0; f < numFields; f++) {
3791 PetscInt fDof;
3792 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr);
3793 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3794 }
3795 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices);CHKERRQ(ierr);
3796 } else {
3797 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices);CHKERRQ(ierr);
3798 rowOffsets[1] = offsetsCopy[0];
3799 }
3800
3801 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr);
3802 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr);
3803 leafEnd = leafStart + numLeaves;
3804 for (l = leafStart; l < leafEnd; l++) {
3805 PetscInt numIndices, childId, offset;
3806 const PetscInt *childIndices;
3807
3808 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr);
3809 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr);
3810 childId = rootIndices[offset++];
3811 childIndices = &rootIndices[offset];
3812 numIndices--;
3813
3814 if (childId == -1) { /* equivalent points: scatter */
3815 PetscInt i;
3816
3817 for (i = 0; i < numIndices; i++) {
3818 PetscInt colIndex = childIndices[i];
3819 PetscInt rowIndex = parentIndices[i];
3820 if (rowIndex < 0) continue;
3821 if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3822 if (colIndex >= colStart && colIndex < colEnd) {
3823 nnzD[rowIndex - rowStart] = 1;
3824 }
3825 else {
3826 nnzO[rowIndex - rowStart] = 1;
3827 }
3828 }
3829 }
3830 else {
3831 PetscInt parentId, f, lim;
3832
3833 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr);
3834
3835 lim = PetscMax(1,numFields);
3836 offsets[0] = 0;
3837 if (numFields) {
3838 PetscInt f;
3839
3840 for (f = 0; f < numFields; f++) {
3841 PetscInt fDof;
3842 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr);
3843
3844 offsets[f + 1] = fDof + offsets[f];
3845 }
3846 }
3847 else {
3848 PetscInt cDof;
3849
3850 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr);
3851 offsets[1] = cDof;
3852 }
3853 for (f = 0; f < lim; f++) {
3854 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3855 PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3856 PetscInt i, numD = 0, numO = 0;
3857
3858 for (i = childStart; i < childEnd; i++) {
3859 PetscInt colIndex = childIndices[i];
3860
3861 if (colIndex < 0) continue;
3862 if (colIndex >= colStart && colIndex < colEnd) {
3863 numD++;
3864 }
3865 else {
3866 numO++;
3867 }
3868 }
3869 for (i = parentStart; i < parentEnd; i++) {
3870 PetscInt rowIndex = parentIndices[i];
3871
3872 if (rowIndex < 0) continue;
3873 nnzD[rowIndex - rowStart] += numD;
3874 nnzO[rowIndex - rowStart] += numO;
3875 }
3876 }
3877 }
3878 }
3879 }
3880 /* preallocate */
3881 ierr = MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);CHKERRQ(ierr);
3882 ierr = PetscFree2(nnzD,nnzO);CHKERRQ(ierr);
3883 /* insert values */
3884 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr);
3885 for (p = pStartC; p < pEndC; p++) {
3886 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3887
3888 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
3889 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr);
3890 if ((dof - cdof) <= 0) continue;
3891 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr);
3892
3893 rowOffsets[0] = 0;
3894 offsetsCopy[0] = 0;
3895 if (numFields) {
3896 PetscInt f;
3897
3898 for (f = 0; f < numFields; f++) {
3899 PetscInt fDof;
3900 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr);
3901 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3902 }
3903 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices);CHKERRQ(ierr);
3904 } else {
3905 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices);CHKERRQ(ierr);
3906 rowOffsets[1] = offsetsCopy[0];
3907 }
3908
3909 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr);
3910 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr);
3911 leafEnd = leafStart + numLeaves;
3912 for (l = leafStart; l < leafEnd; l++) {
3913 PetscInt numIndices, childId, offset;
3914 const PetscInt *childIndices;
3915
3916 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr);
3917 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr);
3918 childId = rootIndices[offset++];
3919 childIndices = &rootIndices[offset];
3920 numIndices--;
3921
3922 if (childId == -1) { /* equivalent points: scatter */
3923 PetscInt i;
3924
3925 for (i = 0; i < numIndices; i++) {
3926 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr);
3927 }
3928 }
3929 else {
3930 PetscInt parentId, f, lim;
3931
3932 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr);
3933
3934 lim = PetscMax(1,numFields);
3935 offsets[0] = 0;
3936 if (numFields) {
3937 PetscInt f;
3938
3939 for (f = 0; f < numFields; f++) {
3940 PetscInt fDof;
3941 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr);
3942
3943 offsets[f + 1] = fDof + offsets[f];
3944 }
3945 }
3946 else {
3947 PetscInt cDof;
3948
3949 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr);
3950 offsets[1] = cDof;
3951 }
3952 for (f = 0; f < lim; f++) {
3953 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3954 PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3955 const PetscInt *colIndices = &childIndices[offsets[f]];
3956
3957 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr);
3958 }
3959 }
3960 }
3961 }
3962 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr);
3963 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr);
3964 ierr = PetscFree(parentIndices);CHKERRQ(ierr);
3965 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr);
3966 ierr = PetscFree(rootIndices);CHKERRQ(ierr);
3967 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr);
3968
3969 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3970 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3971 PetscFunctionReturn(0);
3972 }
3973
DMPlexTransferVecTree_Interpolate(DM coarse,Vec vecCoarseLocal,DM fine,Vec vecFine,PetscSF coarseToFine,PetscInt * cids,Vec grad,Vec cellGeom)3974 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3975 {
3976 PetscSF coarseToFineEmbedded;
3977 PetscSection globalCoarse, globalFine;
3978 PetscSection localCoarse, localFine;
3979 PetscSection aSec, cSec;
3980 PetscSection rootValuesSec;
3981 PetscSection leafValuesSec;
3982 PetscScalar *rootValues, *leafValues;
3983 IS aIS;
3984 const PetscInt *anchors;
3985 Mat cMat;
3986 PetscInt numFields;
3987 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3988 PetscInt aStart, aEnd, cStart, cEnd;
3989 PetscInt *maxChildIds;
3990 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3991 PetscFV fv = NULL;
3992 PetscInt dim, numFVcomps = -1, fvField = -1;
3993 DM cellDM = NULL, gradDM = NULL;
3994 const PetscScalar *cellGeomArray = NULL;
3995 const PetscScalar *gradArray = NULL;
3996 PetscErrorCode ierr;
3997
3998 PetscFunctionBegin;
3999 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr);
4000 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
4001 ierr = DMPlexGetSimplexOrBoxCells(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr);
4002 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr);
4003 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr);
4004 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr);
4005 { /* winnow fine points that don't have global dofs out of the sf */
4006 PetscInt nleaves, l;
4007 const PetscInt *leaves;
4008 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
4009
4010 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr);
4011
4012 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
4013 PetscInt p = leaves ? leaves[l] : l;
4014
4015 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
4016 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
4017 if ((dof - cdof) > 0) {
4018 numPointsWithDofs++;
4019 }
4020 }
4021 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr);
4022 for (l = 0, offset = 0; l < nleaves; l++) {
4023 PetscInt p = leaves ? leaves[l] : l;
4024
4025 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr);
4026 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr);
4027 if ((dof - cdof) > 0) {
4028 pointsWithDofs[offset++] = l;
4029 }
4030 }
4031 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr);
4032 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr);
4033 }
4034 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
4035 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr);
4036 for (p = pStartC; p < pEndC; p++) {
4037 maxChildIds[p - pStartC] = -2;
4038 }
4039 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr);
4040 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr);
4041
4042 ierr = DMGetLocalSection(coarse,&localCoarse);CHKERRQ(ierr);
4043 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr);
4044
4045 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr);
4046 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr);
4047 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
4048
4049 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr);
4050 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr);
4051
4052 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
4053 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr);
4054 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr);
4055 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr);
4056 {
4057 PetscInt maxFields = PetscMax(1,numFields) + 1;
4058 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr);
4059 }
4060 if (grad) {
4061 PetscInt i;
4062
4063 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr);
4064 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr);
4065 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr);
4066 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr);
4067 for (i = 0; i < PetscMax(1,numFields); i++) {
4068 PetscObject obj;
4069 PetscClassId id;
4070
4071 ierr = DMGetField(coarse, i, NULL, &obj);CHKERRQ(ierr);
4072 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr);
4073 if (id == PETSCFV_CLASSID) {
4074 fv = (PetscFV) obj;
4075 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr);
4076 fvField = i;
4077 break;
4078 }
4079 }
4080 }
4081
4082 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
4083 PetscInt dof;
4084 PetscInt maxChildId = maxChildIds[p - pStartC];
4085 PetscInt numValues = 0;
4086
4087 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
4088 if (dof < 0) {
4089 dof = -(dof + 1);
4090 }
4091 offsets[0] = 0;
4092 newOffsets[0] = 0;
4093 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
4094 PetscInt *closure = NULL, closureSize, cl;
4095
4096 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
4097 for (cl = 0; cl < closureSize; cl++) { /* get the closure */
4098 PetscInt c = closure[2 * cl], clDof;
4099
4100 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr);
4101 numValues += clDof;
4102 }
4103 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
4104 }
4105 else if (maxChildId == -1) {
4106 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr);
4107 }
4108 /* we will pack the column indices with the field offsets */
4109 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
4110 /* also send the centroid, and the gradient */
4111 numValues += dim * (1 + numFVcomps);
4112 }
4113 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr);
4114 }
4115 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr);
4116 {
4117 PetscInt numRootValues;
4118 const PetscScalar *coarseArray;
4119
4120 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr);
4121 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr);
4122 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr);
4123 for (p = pStartC; p < pEndC; p++) {
4124 PetscInt numValues;
4125 PetscInt pValOff;
4126 PetscScalar *pVal;
4127 PetscInt maxChildId = maxChildIds[p - pStartC];
4128
4129 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr);
4130 if (!numValues) {
4131 continue;
4132 }
4133 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr);
4134 pVal = &(rootValues[pValOff]);
4135 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
4136 PetscInt closureSize = numValues;
4137 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr);
4138 if (grad && p >= cellStart && p < cellEnd) {
4139 PetscFVCellGeom *cg;
4140 PetscScalar *gradVals = NULL;
4141 PetscInt i;
4142
4143 pVal += (numValues - dim * (1 + numFVcomps));
4144
4145 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr);
4146 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
4147 pVal += dim;
4148 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr);
4149 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
4150 }
4151 }
4152 else if (maxChildId == -1) {
4153 PetscInt lDof, lOff, i;
4154
4155 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr);
4156 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr);
4157 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
4158 }
4159 }
4160 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr);
4161 ierr = PetscFree(maxChildIds);CHKERRQ(ierr);
4162 }
4163 {
4164 PetscSF valuesSF;
4165 PetscInt *remoteOffsetsValues, numLeafValues;
4166
4167 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr);
4168 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr);
4169 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr);
4170 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr);
4171 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr);
4172 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr);
4173 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr);
4174 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr);
4175 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr);
4176 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr);
4177 ierr = PetscFree(rootValues);CHKERRQ(ierr);
4178 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr);
4179 }
4180 ierr = DMGetLocalSection(fine,&localFine);CHKERRQ(ierr);
4181 {
4182 PetscInt maxDof;
4183 PetscInt *rowIndices;
4184 DM refTree;
4185 PetscInt **refPointFieldN;
4186 PetscScalar ***refPointFieldMats;
4187 PetscSection refConSec, refAnSec;
4188 PetscInt pRefStart,pRefEnd,leafStart,leafEnd;
4189 PetscScalar *pointWork;
4190
4191 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr);
4192 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr);
4193 ierr = DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr);
4194 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr);
4195 ierr = DMCopyDisc(fine,refTree);CHKERRQ(ierr);
4196 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
4197 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr);
4198 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr);
4199 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr);
4200 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr);
4201 ierr = DMPlexGetSimplexOrBoxCells(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr);
4202 for (p = leafStart; p < leafEnd; p++) {
4203 PetscInt gDof, gcDof, gOff, lDof;
4204 PetscInt numValues, pValOff;
4205 PetscInt childId;
4206 const PetscScalar *pVal;
4207 const PetscScalar *fvGradData = NULL;
4208
4209 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr);
4210 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr);
4211 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr);
4212 if ((gDof - gcDof) <= 0) {
4213 continue;
4214 }
4215 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr);
4216 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr);
4217 if (!numValues) continue;
4218 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr);
4219 pVal = &leafValues[pValOff];
4220 offsets[0] = 0;
4221 offsetsCopy[0] = 0;
4222 newOffsets[0] = 0;
4223 newOffsetsCopy[0] = 0;
4224 childId = cids[p - pStartF];
4225 if (numFields) {
4226 PetscInt f;
4227 for (f = 0; f < numFields; f++) {
4228 PetscInt rowDof;
4229
4230 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr);
4231 offsets[f + 1] = offsets[f] + rowDof;
4232 offsetsCopy[f + 1] = offsets[f + 1];
4233 /* TODO: closure indices */
4234 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4235 }
4236 ierr = DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,rowIndices);CHKERRQ(ierr);
4237 }
4238 else {
4239 offsets[0] = 0;
4240 offsets[1] = lDof;
4241 newOffsets[0] = 0;
4242 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4243 ierr = DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,rowIndices);CHKERRQ(ierr);
4244 }
4245 if (childId == -1) { /* no child interpolation: one nnz per */
4246 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr);
4247 } else {
4248 PetscInt f;
4249
4250 if (grad && p >= cellStart && p < cellEnd) {
4251 numValues -= (dim * (1 + numFVcomps));
4252 fvGradData = &pVal[numValues];
4253 }
4254 for (f = 0; f < PetscMax(1,numFields); f++) {
4255 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4256 PetscInt numRows = offsets[f+1] - offsets[f];
4257 PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4258 const PetscScalar *cVal = &pVal[newOffsets[f]];
4259 PetscScalar *rVal = &pointWork[offsets[f]];
4260 PetscInt i, j;
4261
4262 #if 0
4263 ierr = PetscInfo5(coarse,"childId %D, numRows %D, numCols %D, refPointFieldN %D maxDof %D\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof);CHKERRQ(ierr);
4264 #endif
4265 for (i = 0; i < numRows; i++) {
4266 PetscScalar val = 0.;
4267 for (j = 0; j < numCols; j++) {
4268 val += childMat[i * numCols + j] * cVal[j];
4269 }
4270 rVal[i] = val;
4271 }
4272 if (f == fvField && p >= cellStart && p < cellEnd) {
4273 PetscReal centroid[3];
4274 PetscScalar diff[3];
4275 const PetscScalar *parentCentroid = &fvGradData[0];
4276 const PetscScalar *gradient = &fvGradData[dim];
4277
4278 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr);
4279 for (i = 0; i < dim; i++) {
4280 diff[i] = centroid[i] - parentCentroid[i];
4281 }
4282 for (i = 0; i < numFVcomps; i++) {
4283 PetscScalar val = 0.;
4284
4285 for (j = 0; j < dim; j++) {
4286 val += gradient[dim * i + j] * diff[j];
4287 }
4288 rVal[i] += val;
4289 }
4290 }
4291 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr);
4292 }
4293 }
4294 }
4295 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr);
4296 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr);
4297 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr);
4298 }
4299 ierr = PetscFree(leafValues);CHKERRQ(ierr);
4300 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr);
4301 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr);
4302 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr);
4303 PetscFunctionReturn(0);
4304 }
4305
DMPlexTransferVecTree_Inject(DM fine,Vec vecFine,DM coarse,Vec vecCoarse,PetscSF coarseToFine,PetscInt * cids)4306 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4307 {
4308 DM refTree;
4309 PetscSection multiRootSec, rootIndicesSec;
4310 PetscSection globalCoarse, globalFine;
4311 PetscSection localCoarse, localFine;
4312 PetscSection cSecRef;
4313 PetscInt *parentIndices, pRefStart, pRefEnd;
4314 PetscScalar *rootValues, *parentValues;
4315 Mat injRef;
4316 PetscInt numFields, maxDof;
4317 PetscInt pStartC, pEndC, pStartF, pEndF, p;
4318 PetscInt *offsets, *offsetsCopy, *rowOffsets;
4319 PetscLayout rowMap, colMap;
4320 PetscInt rowStart, rowEnd, colStart, colEnd;
4321 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4322 PetscErrorCode ierr;
4323
4324 PetscFunctionBegin;
4325
4326 /* get the templates for the fine-to-coarse injection from the reference tree */
4327 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr);
4328 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr);
4329 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr);
4330 ierr = DMCopyDisc(coarse,refTree);CHKERRQ(ierr);
4331 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr);
4332 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr);
4333 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr);
4334
4335 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr);
4336 ierr = DMGetLocalSection(fine,&localFine);CHKERRQ(ierr);
4337 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr);
4338 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr);
4339 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr);
4340 ierr = DMGetLocalSection(coarse,&localCoarse);CHKERRQ(ierr);
4341 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr);
4342 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr);
4343 {
4344 PetscInt maxFields = PetscMax(1,numFields) + 1;
4345 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr);
4346 }
4347
4348 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr);
4349
4350 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr);
4351
4352 /* count indices */
4353 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr);
4354 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr);
4355 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr);
4356 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr);
4357 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr);
4358 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr);
4359 /* insert values */
4360 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr);
4361 for (p = pStartC; p < pEndC; p++) {
4362 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4363 PetscBool contribute = PETSC_FALSE;
4364
4365 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr);
4366 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr);
4367 if ((dof - cdof) <= 0) continue;
4368 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr);
4369 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr);
4370
4371 rowOffsets[0] = 0;
4372 offsetsCopy[0] = 0;
4373 if (numFields) {
4374 PetscInt f;
4375
4376 for (f = 0; f < numFields; f++) {
4377 PetscInt fDof;
4378 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr);
4379 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4380 }
4381 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,parentIndices);CHKERRQ(ierr);
4382 } else {
4383 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,parentIndices);CHKERRQ(ierr);
4384 rowOffsets[1] = offsetsCopy[0];
4385 }
4386
4387 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr);
4388 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr);
4389 leafEnd = leafStart + numLeaves;
4390 for (l = 0; l < dof; l++) parentValues[l] = 0.;
4391 for (l = leafStart; l < leafEnd; l++) {
4392 PetscInt numIndices, childId, offset;
4393 const PetscScalar *childValues;
4394
4395 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr);
4396 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr);
4397 childId = (PetscInt) PetscRealPart(rootValues[offset++]);
4398 childValues = &rootValues[offset];
4399 numIndices--;
4400
4401 if (childId == -2) { /* skip */
4402 continue;
4403 } else if (childId == -1) { /* equivalent points: scatter */
4404 PetscInt m;
4405
4406 contribute = PETSC_TRUE;
4407 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4408 } else { /* contributions from children: sum with injectors from reference tree */
4409 PetscInt parentId, f, lim;
4410
4411 contribute = PETSC_TRUE;
4412 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr);
4413
4414 lim = PetscMax(1,numFields);
4415 offsets[0] = 0;
4416 if (numFields) {
4417 PetscInt f;
4418
4419 for (f = 0; f < numFields; f++) {
4420 PetscInt fDof;
4421 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr);
4422
4423 offsets[f + 1] = fDof + offsets[f];
4424 }
4425 }
4426 else {
4427 PetscInt cDof;
4428
4429 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr);
4430 offsets[1] = cDof;
4431 }
4432 for (f = 0; f < lim; f++) {
4433 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4434 PetscInt n = offsets[f+1]-offsets[f];
4435 PetscInt m = rowOffsets[f+1]-rowOffsets[f];
4436 PetscInt i, j;
4437 const PetscScalar *colValues = &childValues[offsets[f]];
4438
4439 for (i = 0; i < m; i++) {
4440 PetscScalar val = 0.;
4441 for (j = 0; j < n; j++) {
4442 val += childMat[n * i + j] * colValues[j];
4443 }
4444 parentValues[rowOffsets[f] + i] += val;
4445 }
4446 }
4447 }
4448 }
4449 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);}
4450 }
4451 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr);
4452 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr);
4453 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr);
4454 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr);
4455 ierr = PetscFree(rootValues);CHKERRQ(ierr);
4456 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr);
4457 PetscFunctionReturn(0);
4458 }
4459
4460 /*@
4461 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4462 that can be represented by a common reference tree used by both. This routine can be used for a combination of
4463 coarsening and refinement at the same time.
4464
4465 collective
4466
4467 Input Parameters:
4468 + dmIn - The DMPlex mesh for the input vector
4469 . vecIn - The input vector
4470 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4471 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4472 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4473 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4474 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies
4475 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4476 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly
4477 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this
4478 point j in dmOut is not a leaf of sfRefine.
4479 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies
4480 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4481 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4482 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4483 - time - Used if boundary values are time dependent.
4484
4485 Output Parameters:
4486 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
4487 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume
4488 method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4489 coarse points to fine points.
4490
4491 Level: developer
4492
4493 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients()
4494 @*/
DMPlexTransferVecTree(DM dmIn,Vec vecIn,DM dmOut,Vec vecOut,PetscSF sfRefine,PetscSF sfCoarsen,PetscInt * cidsRefine,PetscInt * cidsCoarsen,PetscBool useBCs,PetscReal time)4495 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4496 {
4497 PetscErrorCode ierr;
4498
4499 PetscFunctionBegin;
4500 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr);
4501 if (sfRefine) {
4502 Vec vecInLocal;
4503 DM dmGrad = NULL;
4504 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4505
4506 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr);
4507 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr);
4508 {
4509 PetscInt numFields, i;
4510
4511 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr);
4512 for (i = 0; i < numFields; i++) {
4513 PetscObject obj;
4514 PetscClassId classid;
4515
4516 ierr = DMGetField(dmIn, i, NULL, &obj);CHKERRQ(ierr);
4517 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr);
4518 if (classid == PETSCFV_CLASSID) {
4519 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr);
4520 break;
4521 }
4522 }
4523 }
4524 if (useBCs) {
4525 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr);
4526 }
4527 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr);
4528 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr);
4529 if (dmGrad) {
4530 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr);
4531 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr);
4532 }
4533 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr);
4534 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr);
4535 if (dmGrad) {
4536 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr);
4537 }
4538 }
4539 if (sfCoarsen) {
4540 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr);
4541 }
4542 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr);
4543 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr);
4544 PetscFunctionReturn(0);
4545 }
4546