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