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); 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 ierr = DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 357 ierr = DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 358 ierr = DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 359 ierr = 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 ierr = 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 ierr = PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 414 ierr = 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 = PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE);CHKERRQ(ierr); 1266 ierr = PetscSpaceSetNumComponents(bspace,Nc);CHKERRQ(ierr); 1267 ierr = PetscSpaceSetNumVariables(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 const PetscReal xi0[3] = {-1.,-1.,-1.}; 1340 1341 CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp); 1342 CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]); 1343 } 1344 ierr = EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);CHKERRQ(ierr); 1345 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1346 ierr = MatDenseGetArray(Xmat,&X);CHKERRQ(ierr); 1347 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1348 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1349 childOffsets[0] = 0; 1350 for (i = 0; i < closureSize; i++) { 1351 PetscInt p = closure[2*i]; 1352 PetscInt dof; 1353 1354 if (numFields) { 1355 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1356 } 1357 else { 1358 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1359 } 1360 childOffsets[i+1]=childOffsets[i]+dof; 1361 } 1362 parentOffsets[0] = 0; 1363 for (i = 0; i < closureSizeP; i++) { 1364 PetscInt p = closureP[2*i]; 1365 PetscInt dof; 1366 1367 if (numFields) { 1368 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1369 } 1370 else { 1371 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1372 } 1373 parentOffsets[i+1]=parentOffsets[i]+dof; 1374 } 1375 for (i = 0; i < closureSize; i++) { 1376 PetscInt conDof, conOff, aDof, aOff, nWork; 1377 PetscInt p = closure[2*i]; 1378 PetscInt o = closure[2*i+1]; 1379 const PetscInt *perm; 1380 const PetscScalar *flip; 1381 1382 if (p < conStart || p >= conEnd) continue; 1383 if (numFields) { 1384 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1385 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1386 } 1387 else { 1388 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1389 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1390 } 1391 if (!conDof) continue; 1392 perm = (perms && perms[i]) ? perms[i][o] : NULL; 1393 flip = (flips && flips[i]) ? flips[i][o] : NULL; 1394 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1395 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1396 nWork = childOffsets[i+1]-childOffsets[i]; 1397 for (k = 0; k < aDof; k++) { 1398 PetscInt a = anchors[aOff + k]; 1399 PetscInt aSecDof, aSecOff; 1400 1401 if (numFields) { 1402 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1403 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1404 } 1405 else { 1406 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1407 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1408 } 1409 if (!aSecDof) continue; 1410 1411 for (j = 0; j < closureSizeP; j++) { 1412 PetscInt q = closureP[2*j]; 1413 PetscInt oq = closureP[2*j+1]; 1414 1415 if (q == a) { 1416 PetscInt r, s, nWorkP; 1417 const PetscInt *permP; 1418 const PetscScalar *flipP; 1419 1420 permP = (perms && perms[j]) ? perms[j][oq] : NULL; 1421 flipP = (flips && flips[j]) ? flips[j][oq] : NULL; 1422 nWorkP = parentOffsets[j+1]-parentOffsets[j]; 1423 /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the 1424 * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArray is 1425 * column-major, so transpose-transpose = do nothing */ 1426 for (r = 0; r < nWork; r++) { 1427 for (s = 0; s < nWorkP; s++) { 1428 scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])]; 1429 } 1430 } 1431 for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + r;} 1432 for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;} 1433 if (flip) { 1434 for (r = 0; r < nWork; r++) { 1435 for (s = 0; s < nWorkP; s++) { 1436 scwork[r * nWorkP + s] *= flip[r]; 1437 } 1438 } 1439 } 1440 if (flipP) { 1441 for (r = 0; r < nWork; r++) { 1442 for (s = 0; s < nWorkP; s++) { 1443 scwork[r * nWorkP + s] *= flipP[s]; 1444 } 1445 } 1446 } 1447 ierr = MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);CHKERRQ(ierr); 1448 break; 1449 } 1450 } 1451 } 1452 } 1453 ierr = MatDenseRestoreArray(Xmat,&X);CHKERRQ(ierr); 1454 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1455 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1456 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1457 } 1458 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1459 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1460 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1461 ierr = PetscFree(scwork);CHKERRQ(ierr); 1462 ierr = PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);CHKERRQ(ierr); 1463 if (id == PETSCFV_CLASSID) { 1464 ierr = PetscSpaceDestroy(&bspace);CHKERRQ(ierr); 1465 } 1466 } 1467 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1468 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1469 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1470 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1471 1472 PetscFunctionReturn(0); 1473 } 1474 1475 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1476 { 1477 Mat refCmat; 1478 PetscDS ds; 1479 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN; 1480 PetscScalar ***refPointFieldMats; 1481 PetscSection refConSec, refAnSec, refSection; 1482 IS refAnIS; 1483 const PetscInt *refAnchors; 1484 const PetscInt **perms; 1485 const PetscScalar **flips; 1486 PetscErrorCode ierr; 1487 1488 PetscFunctionBegin; 1489 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1490 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1491 maxFields = PetscMax(1,numFields); 1492 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1493 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1494 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1495 ierr = DMGetSection(refTree,&refSection);CHKERRQ(ierr); 1496 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1497 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 1498 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr); 1499 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1500 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1501 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 1502 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr); 1503 for (p = pRefStart; p < pRefEnd; p++) { 1504 PetscInt parent, closureSize, *closure = NULL, pDof; 1505 1506 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1507 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1508 if (!pDof || parent == p) continue; 1509 1510 ierr = PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 1511 ierr = PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr); 1512 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1513 for (f = 0; f < maxFields; f++) { 1514 PetscInt cDof, cOff, numCols, r, i; 1515 1516 if (f < numFields) { 1517 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1518 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 1519 ierr = PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1520 } else { 1521 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1522 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 1523 ierr = PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1524 } 1525 1526 for (r = 0; r < cDof; r++) { 1527 rows[r] = cOff + r; 1528 } 1529 numCols = 0; 1530 for (i = 0; i < closureSize; i++) { 1531 PetscInt q = closure[2*i]; 1532 PetscInt aDof, aOff, j; 1533 const PetscInt *perm = perms ? perms[i] : NULL; 1534 1535 if (numFields) { 1536 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1537 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1538 } 1539 else { 1540 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1541 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1542 } 1543 1544 for (j = 0; j < aDof; j++) { 1545 cols[numCols++] = aOff + (perm ? perm[j] : j); 1546 } 1547 } 1548 refPointFieldN[p-pRefStart][f] = numCols; 1549 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1550 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1551 if (flips) { 1552 PetscInt colOff = 0; 1553 1554 for (i = 0; i < closureSize; i++) { 1555 PetscInt q = closure[2*i]; 1556 PetscInt aDof, aOff, j; 1557 const PetscScalar *flip = flips ? flips[i] : NULL; 1558 1559 if (numFields) { 1560 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1561 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1562 } 1563 else { 1564 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1565 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1566 } 1567 if (flip) { 1568 PetscInt k; 1569 for (k = 0; k < cDof; k++) { 1570 for (j = 0; j < aDof; j++) { 1571 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j]; 1572 } 1573 } 1574 } 1575 colOff += aDof; 1576 } 1577 } 1578 if (numFields) { 1579 ierr = PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1580 } else { 1581 ierr = PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1582 } 1583 } 1584 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1585 } 1586 *childrenMats = refPointFieldMats; 1587 *childrenN = refPointFieldN; 1588 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1589 ierr = PetscFree(rows);CHKERRQ(ierr); 1590 ierr = PetscFree(cols);CHKERRQ(ierr); 1591 PetscFunctionReturn(0); 1592 } 1593 1594 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1595 { 1596 PetscDS ds; 1597 PetscInt **refPointFieldN; 1598 PetscScalar ***refPointFieldMats; 1599 PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f; 1600 PetscSection refConSec; 1601 PetscErrorCode ierr; 1602 1603 PetscFunctionBegin; 1604 refPointFieldN = *childrenN; 1605 *childrenN = NULL; 1606 refPointFieldMats = *childrenMats; 1607 *childrenMats = NULL; 1608 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1609 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1610 maxFields = PetscMax(1,numFields); 1611 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 1612 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1613 for (p = pRefStart; p < pRefEnd; p++) { 1614 PetscInt parent, pDof; 1615 1616 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1617 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1618 if (!pDof || parent == p) continue; 1619 1620 for (f = 0; f < maxFields; f++) { 1621 PetscInt cDof; 1622 1623 if (numFields) { 1624 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1625 } 1626 else { 1627 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1628 } 1629 1630 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 1631 } 1632 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 1633 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr); 1634 } 1635 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 1636 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr); 1637 PetscFunctionReturn(0); 1638 } 1639 1640 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat) 1641 { 1642 DM refTree; 1643 PetscDS ds; 1644 Mat refCmat; 1645 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1646 PetscScalar ***refPointFieldMats, *pointWork; 1647 PetscSection refConSec, refAnSec, anSec; 1648 IS refAnIS, anIS; 1649 const PetscInt *anchors; 1650 PetscErrorCode ierr; 1651 1652 PetscFunctionBegin; 1653 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1654 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1655 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1656 maxFields = PetscMax(1,numFields); 1657 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 1658 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 1659 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1660 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1661 ierr = DMPlexGetAnchors(dm,&anSec,&anIS);CHKERRQ(ierr); 1662 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr); 1663 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1664 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr); 1665 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1666 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1667 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr); 1668 1669 /* step 1: get submats for every constrained point in the reference tree */ 1670 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1671 1672 /* step 2: compute the preorder */ 1673 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1674 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr); 1675 for (p = pStart; p < pEnd; p++) { 1676 perm[p - pStart] = p; 1677 iperm[p - pStart] = p-pStart; 1678 } 1679 for (p = 0; p < pEnd - pStart;) { 1680 PetscInt point = perm[p]; 1681 PetscInt parent; 1682 1683 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr); 1684 if (parent == point) { 1685 p++; 1686 } 1687 else { 1688 PetscInt size, closureSize, *closure = NULL, i; 1689 1690 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1691 for (i = 0; i < closureSize; i++) { 1692 PetscInt q = closure[2*i]; 1693 if (iperm[q-pStart] > iperm[point-pStart]) { 1694 /* swap */ 1695 perm[p] = q; 1696 perm[iperm[q-pStart]] = point; 1697 iperm[point-pStart] = iperm[q-pStart]; 1698 iperm[q-pStart] = p; 1699 break; 1700 } 1701 } 1702 size = closureSize; 1703 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1704 if (i == size) { 1705 p++; 1706 } 1707 } 1708 } 1709 1710 /* step 3: fill the constraint matrix */ 1711 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1712 * allow progressive fill without assembly, so we are going to set up the 1713 * values outside of the Mat first. 1714 */ 1715 { 1716 PetscInt nRows, row, nnz; 1717 PetscBool done; 1718 const PetscInt *ia, *ja; 1719 PetscScalar *vals; 1720 1721 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1722 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1723 nnz = ia[nRows]; 1724 /* malloc and then zero rows right before we fill them: this way valgrind 1725 * can tell if we are doing progressive fill in the wrong order */ 1726 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr); 1727 for (p = 0; p < pEnd - pStart; p++) { 1728 PetscInt parent, childid, closureSize, *closure = NULL; 1729 PetscInt point = perm[p], pointDof; 1730 1731 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr); 1732 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1733 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr); 1734 if (!pointDof) continue; 1735 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1736 for (f = 0; f < maxFields; f++) { 1737 PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset; 1738 PetscScalar *pointMat; 1739 const PetscInt **perms; 1740 const PetscScalar **flips; 1741 1742 if (numFields) { 1743 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr); 1744 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr); 1745 } 1746 else { 1747 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr); 1748 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr); 1749 } 1750 if (!cDof) continue; 1751 if (numFields) {ierr = PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);} 1752 else {ierr = PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr);} 1753 1754 /* make sure that every row for this point is the same size */ 1755 #if defined(PETSC_USE_DEBUG) 1756 for (r = 0; r < cDof; r++) { 1757 if (cDof > 1 && r) { 1758 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])); 1759 } 1760 } 1761 #endif 1762 /* zero rows */ 1763 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1764 vals[i] = 0.; 1765 } 1766 matOffset = ia[cOff]; 1767 numFillCols = ia[cOff+1] - matOffset; 1768 pointMat = refPointFieldMats[childid-pRefStart][f]; 1769 numCols = refPointFieldN[childid-pRefStart][f]; 1770 offset = 0; 1771 for (i = 0; i < closureSize; i++) { 1772 PetscInt q = closure[2*i]; 1773 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1774 const PetscInt *perm = perms ? perms[i] : NULL; 1775 const PetscScalar *flip = flips ? flips[i] : NULL; 1776 1777 qConDof = qConOff = 0; 1778 if (numFields) { 1779 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr); 1780 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr); 1781 if (q >= conStart && q < conEnd) { 1782 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr); 1783 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr); 1784 } 1785 } 1786 else { 1787 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr); 1788 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr); 1789 if (q >= conStart && q < conEnd) { 1790 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr); 1791 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr); 1792 } 1793 } 1794 if (!aDof) continue; 1795 if (qConDof) { 1796 /* this point has anchors: its rows of the matrix should already 1797 * be filled, thanks to preordering */ 1798 /* first multiply into pointWork, then set in matrix */ 1799 PetscInt aMatOffset = ia[qConOff]; 1800 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1801 for (r = 0; r < cDof; r++) { 1802 for (j = 0; j < aNumFillCols; j++) { 1803 PetscScalar inVal = 0; 1804 for (k = 0; k < aDof; k++) { 1805 PetscInt col = perm ? perm[k] : k; 1806 1807 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.); 1808 } 1809 pointWork[r * aNumFillCols + j] = inVal; 1810 } 1811 } 1812 /* assume that the columns are sorted, spend less time searching */ 1813 for (j = 0, k = 0; j < aNumFillCols; j++) { 1814 PetscInt col = ja[aMatOffset + j]; 1815 for (;k < numFillCols; k++) { 1816 if (ja[matOffset + k] == col) { 1817 break; 1818 } 1819 } 1820 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col); 1821 for (r = 0; r < cDof; r++) { 1822 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1823 } 1824 } 1825 } 1826 else { 1827 /* find where to put this portion of pointMat into the matrix */ 1828 for (k = 0; k < numFillCols; k++) { 1829 if (ja[matOffset + k] == aOff) { 1830 break; 1831 } 1832 } 1833 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff); 1834 for (r = 0; r < cDof; r++) { 1835 for (j = 0; j < aDof; j++) { 1836 PetscInt col = perm ? perm[j] : j; 1837 1838 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.); 1839 } 1840 } 1841 } 1842 offset += aDof; 1843 } 1844 if (numFields) { 1845 ierr = PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1846 } else { 1847 ierr = PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1848 } 1849 } 1850 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1851 } 1852 for (row = 0; row < nRows; row++) { 1853 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr); 1854 } 1855 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1856 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1857 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1858 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1859 ierr = PetscFree(vals);CHKERRQ(ierr); 1860 } 1861 1862 /* clean up */ 1863 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr); 1864 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr); 1865 ierr = PetscFree(pointWork);CHKERRQ(ierr); 1866 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1867 PetscFunctionReturn(0); 1868 } 1869 1870 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of 1871 * a non-conforming mesh. Local refinement comes later */ 1872 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm) 1873 { 1874 DM K; 1875 PetscMPIInt rank; 1876 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd; 1877 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations; 1878 PetscInt *Kembedding; 1879 PetscInt *cellClosure=NULL, nc; 1880 PetscScalar *newVertexCoords; 1881 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset; 1882 PetscSection parentSection; 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1887 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1888 ierr = DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);CHKERRQ(ierr); 1889 ierr = DMSetDimension(*ncdm,dim);CHKERRQ(ierr); 1890 1891 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1892 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);CHKERRQ(ierr); 1893 ierr = DMPlexGetReferenceTree(dm,&K);CHKERRQ(ierr); 1894 if (!rank) { 1895 /* compute the new charts */ 1896 ierr = PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);CHKERRQ(ierr); 1897 offset = 0; 1898 for (d = 0; d <= dim; d++) { 1899 PetscInt pOldCount, kStart, kEnd, k; 1900 1901 pNewStart[d] = offset; 1902 ierr = DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);CHKERRQ(ierr); 1903 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1904 pOldCount = pOldEnd[d] - pOldStart[d]; 1905 /* adding the new points */ 1906 pNewCount[d] = pOldCount + kEnd - kStart; 1907 if (!d) { 1908 /* removing the cell */ 1909 pNewCount[d]--; 1910 } 1911 for (k = kStart; k < kEnd; k++) { 1912 PetscInt parent; 1913 ierr = DMPlexGetTreeParent(K,k,&parent,NULL);CHKERRQ(ierr); 1914 if (parent == k) { 1915 /* avoid double counting points that won't actually be new */ 1916 pNewCount[d]--; 1917 } 1918 } 1919 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1920 offset = pNewEnd[d]; 1921 1922 } 1923 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]); 1924 /* get the current closure of the cell that we are removing */ 1925 ierr = DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 1926 1927 ierr = PetscMalloc1(pNewEnd[dim],&newConeSizes);CHKERRQ(ierr); 1928 { 1929 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1930 1931 ierr = DMPlexGetChart(K,&kStart,&kEnd);CHKERRQ(ierr); 1932 ierr = PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);CHKERRQ(ierr); 1933 1934 for (k = kStart; k < kEnd; k++) { 1935 perm[k - kStart] = k; 1936 iperm [k - kStart] = k - kStart; 1937 preOrient[k - kStart] = 0; 1938 } 1939 1940 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1941 for (j = 1; j < closureSizeK; j++) { 1942 PetscInt parentOrientA = closureK[2*j+1]; 1943 PetscInt parentOrientB = cellClosure[2*j+1]; 1944 PetscInt p, q; 1945 1946 p = closureK[2*j]; 1947 q = cellClosure[2*j]; 1948 for (d = 0; d <= dim; d++) { 1949 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1950 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1951 } 1952 } 1953 if (parentOrientA != parentOrientB) { 1954 PetscInt numChildren, i; 1955 const PetscInt *children; 1956 1957 ierr = DMPlexGetTreeChildren(K,p,&numChildren,&children);CHKERRQ(ierr); 1958 for (i = 0; i < numChildren; i++) { 1959 PetscInt kPerm, oPerm; 1960 1961 k = children[i]; 1962 ierr = DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);CHKERRQ(ierr); 1963 /* perm = what refTree position I'm in */ 1964 perm[kPerm-kStart] = k; 1965 /* iperm = who is at this position */ 1966 iperm[k-kStart] = kPerm-kStart; 1967 preOrient[kPerm-kStart] = oPerm; 1968 } 1969 } 1970 } 1971 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1972 } 1973 ierr = PetscSectionSetChart(parentSection,0,pNewEnd[dim]);CHKERRQ(ierr); 1974 offset = 0; 1975 numNewCones = 0; 1976 for (d = 0; d <= dim; d++) { 1977 PetscInt kStart, kEnd, k; 1978 PetscInt p; 1979 PetscInt size; 1980 1981 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 1982 /* skip cell 0 */ 1983 if (p == cell) continue; 1984 /* old cones to new cones */ 1985 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 1986 newConeSizes[offset++] = size; 1987 numNewCones += size; 1988 } 1989 1990 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1991 for (k = kStart; k < kEnd; k++) { 1992 PetscInt kParent; 1993 1994 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 1995 if (kParent != k) { 1996 Kembedding[k] = offset; 1997 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 1998 newConeSizes[offset++] = size; 1999 numNewCones += size; 2000 if (kParent != 0) { 2001 ierr = PetscSectionSetDof(parentSection,Kembedding[k],1);CHKERRQ(ierr); 2002 } 2003 } 2004 } 2005 } 2006 2007 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2008 ierr = PetscSectionGetStorageSize(parentSection,&numPointsWithParents);CHKERRQ(ierr); 2009 ierr = PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);CHKERRQ(ierr); 2010 ierr = PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);CHKERRQ(ierr); 2011 2012 /* fill new cones */ 2013 offset = 0; 2014 for (d = 0; d <= dim; d++) { 2015 PetscInt kStart, kEnd, k, l; 2016 PetscInt p; 2017 PetscInt size; 2018 const PetscInt *cone, *orientation; 2019 2020 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2021 /* skip cell 0 */ 2022 if (p == cell) continue; 2023 /* old cones to new cones */ 2024 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 2025 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 2026 ierr = DMPlexGetConeOrientation(dm,p,&orientation);CHKERRQ(ierr); 2027 for (l = 0; l < size; l++) { 2028 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 2029 newOrientations[offset++] = orientation[l]; 2030 } 2031 } 2032 2033 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 2034 for (k = kStart; k < kEnd; k++) { 2035 PetscInt kPerm = perm[k], kParent; 2036 PetscInt preO = preOrient[k]; 2037 2038 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 2039 if (kParent != k) { 2040 /* embed new cones */ 2041 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 2042 ierr = DMPlexGetCone(K,kPerm,&cone);CHKERRQ(ierr); 2043 ierr = DMPlexGetConeOrientation(K,kPerm,&orientation);CHKERRQ(ierr); 2044 for (l = 0; l < size; l++) { 2045 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 2046 PetscInt newO, lSize, oTrue; 2047 2048 q = iperm[cone[m]]; 2049 newCones[offset] = Kembedding[q]; 2050 ierr = DMPlexGetConeSize(K,q,&lSize);CHKERRQ(ierr); 2051 oTrue = orientation[m]; 2052 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 2053 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 2054 newOrientations[offset++] = newO; 2055 } 2056 if (kParent != 0) { 2057 PetscInt newPoint = Kembedding[kParent]; 2058 ierr = PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);CHKERRQ(ierr); 2059 parents[pOffset] = newPoint; 2060 childIDs[pOffset] = k; 2061 } 2062 } 2063 } 2064 } 2065 2066 ierr = PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);CHKERRQ(ierr); 2067 2068 /* fill coordinates */ 2069 offset = 0; 2070 { 2071 PetscInt kStart, kEnd, l; 2072 PetscSection vSection; 2073 PetscInt v; 2074 Vec coords; 2075 PetscScalar *coordvals; 2076 PetscInt dof, off; 2077 PetscReal v0[3], J[9], detJ; 2078 2079 #if defined(PETSC_USE_DEBUG) 2080 { 2081 PetscInt k; 2082 ierr = DMPlexGetHeightStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2083 for (k = kStart; k < kEnd; k++) { 2084 ierr = DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2085 if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k); 2086 } 2087 } 2088 #endif 2089 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2090 ierr = DMGetCoordinateSection(dm,&vSection);CHKERRQ(ierr); 2091 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2092 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2093 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 2094 2095 ierr = PetscSectionGetDof(vSection,v,&dof);CHKERRQ(ierr); 2096 ierr = PetscSectionGetOffset(vSection,v,&off);CHKERRQ(ierr); 2097 for (l = 0; l < dof; l++) { 2098 newVertexCoords[offset++] = coordvals[off + l]; 2099 } 2100 } 2101 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2102 2103 ierr = DMGetCoordinateSection(K,&vSection);CHKERRQ(ierr); 2104 ierr = DMGetCoordinatesLocal(K,&coords);CHKERRQ(ierr); 2105 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2106 ierr = DMPlexGetDepthStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2107 for (v = kStart; v < kEnd; v++) { 2108 PetscReal coord[3], newCoord[3]; 2109 PetscInt vPerm = perm[v]; 2110 PetscInt kParent; 2111 const PetscReal xi0[3] = {-1.,-1.,-1.}; 2112 2113 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr); 2114 if (kParent != v) { 2115 /* this is a new vertex */ 2116 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr); 2117 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2118 CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord); 2119 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2120 offset += dim; 2121 } 2122 } 2123 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2124 } 2125 2126 /* need to reverse the order of pNewCount: vertices first, cells last */ 2127 for (d = 0; d < (dim + 1) / 2; d++) { 2128 PetscInt tmp; 2129 2130 tmp = pNewCount[d]; 2131 pNewCount[d] = pNewCount[dim - d]; 2132 pNewCount[dim - d] = tmp; 2133 } 2134 2135 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 2136 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2137 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 2138 2139 /* clean up */ 2140 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 2141 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 2142 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 2143 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 2144 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 2145 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 2146 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 2147 } 2148 else { 2149 PetscInt p, counts[4]; 2150 PetscInt *coneSizes, *cones, *orientations; 2151 Vec coordVec; 2152 PetscScalar *coords; 2153 2154 for (d = 0; d <= dim; d++) { 2155 PetscInt dStart, dEnd; 2156 2157 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 2158 counts[d] = dEnd - dStart; 2159 } 2160 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 2161 for (p = pStart; p < pEnd; p++) { 2162 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 2163 } 2164 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2165 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2166 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2167 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2168 2169 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2170 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2171 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2172 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2173 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2174 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2175 } 2176 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2177 2178 PetscFunctionReturn(0); 2179 } 2180 2181 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2182 { 2183 PetscSF coarseToFineEmbedded; 2184 PetscSection globalCoarse, globalFine; 2185 PetscSection localCoarse, localFine; 2186 PetscSection aSec, cSec; 2187 PetscSection rootIndicesSec, rootMatricesSec; 2188 PetscSection leafIndicesSec, leafMatricesSec; 2189 PetscInt *rootIndices, *leafIndices; 2190 PetscScalar *rootMatrices, *leafMatrices; 2191 IS aIS; 2192 const PetscInt *anchors; 2193 Mat cMat; 2194 PetscInt numFields, maxFields; 2195 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2196 PetscInt aStart, aEnd, cStart, cEnd; 2197 PetscInt *maxChildIds; 2198 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2199 const PetscInt ***perms; 2200 const PetscScalar ***flips; 2201 PetscErrorCode ierr; 2202 2203 PetscFunctionBegin; 2204 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 2205 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 2206 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr); 2207 { /* winnow fine points that don't have global dofs out of the sf */ 2208 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2209 const PetscInt *leaves; 2210 2211 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 2212 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2213 p = leaves ? leaves[l] : l; 2214 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2215 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2216 if ((dof - cdof) > 0) { 2217 numPointsWithDofs++; 2218 } 2219 } 2220 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 2221 for (l = 0, offset = 0; l < nleaves; l++) { 2222 p = leaves ? leaves[l] : l; 2223 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2224 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2225 if ((dof - cdof) > 0) { 2226 pointsWithDofs[offset++] = l; 2227 } 2228 } 2229 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 2230 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 2231 } 2232 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2233 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 2234 for (p = pStartC; p < pEndC; p++) { 2235 maxChildIds[p - pStartC] = -2; 2236 } 2237 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2238 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2239 2240 ierr = DMGetSection(coarse,&localCoarse);CHKERRQ(ierr); 2241 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 2242 2243 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 2244 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 2245 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 2246 2247 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 2248 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 2249 2250 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2251 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 2252 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr); 2253 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr); 2254 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr); 2255 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 2256 maxFields = PetscMax(1,numFields); 2257 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); 2258 ierr = PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);CHKERRQ(ierr); 2259 ierr = PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));CHKERRQ(ierr); 2260 ierr = PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));CHKERRQ(ierr); 2261 2262 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2263 PetscInt dof, matSize = 0; 2264 PetscInt aDof = 0; 2265 PetscInt cDof = 0; 2266 PetscInt maxChildId = maxChildIds[p - pStartC]; 2267 PetscInt numRowIndices = 0; 2268 PetscInt numColIndices = 0; 2269 PetscInt f; 2270 2271 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2272 if (dof < 0) { 2273 dof = -(dof + 1); 2274 } 2275 if (p >= aStart && p < aEnd) { 2276 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2277 } 2278 if (p >= cStart && p < cEnd) { 2279 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr); 2280 } 2281 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2282 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2283 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2284 PetscInt *closure = NULL, closureSize, cl; 2285 2286 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2287 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2288 PetscInt c = closure[2 * cl], clDof; 2289 2290 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 2291 numRowIndices += clDof; 2292 for (f = 0; f < numFields; f++) { 2293 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr); 2294 offsets[f + 1] += clDof; 2295 } 2296 } 2297 for (f = 0; f < numFields; f++) { 2298 offsets[f + 1] += offsets[f]; 2299 newOffsets[f + 1] = offsets[f + 1]; 2300 } 2301 /* get the number of indices needed and their field offsets */ 2302 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2303 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2304 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2305 numColIndices = numRowIndices; 2306 matSize = 0; 2307 } 2308 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2309 matSize = 0; 2310 for (f = 0; f < numFields; f++) { 2311 PetscInt numRow, numCol; 2312 2313 numRow = offsets[f + 1] - offsets[f]; 2314 numCol = newOffsets[f + 1] - newOffsets[f]; 2315 matSize += numRow * numCol; 2316 } 2317 } 2318 else { 2319 matSize = numRowIndices * numColIndices; 2320 } 2321 } else if (maxChildId == -1) { 2322 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2323 PetscInt aOff, a; 2324 2325 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2326 for (f = 0; f < numFields; f++) { 2327 PetscInt fDof; 2328 2329 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2330 offsets[f+1] = fDof; 2331 } 2332 for (a = 0; a < aDof; a++) { 2333 PetscInt anchor = anchors[a + aOff], aLocalDof; 2334 2335 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr); 2336 numColIndices += aLocalDof; 2337 for (f = 0; f < numFields; f++) { 2338 PetscInt fDof; 2339 2340 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2341 newOffsets[f+1] += fDof; 2342 } 2343 } 2344 if (numFields) { 2345 matSize = 0; 2346 for (f = 0; f < numFields; f++) { 2347 matSize += offsets[f+1] * newOffsets[f+1]; 2348 } 2349 } 2350 else { 2351 matSize = numColIndices * dof; 2352 } 2353 } 2354 else { /* no children, and no constraints on dofs: just get the global indices */ 2355 numColIndices = dof; 2356 matSize = 0; 2357 } 2358 } 2359 /* we will pack the column indices with the field offsets */ 2360 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr); 2361 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr); 2362 } 2363 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr); 2364 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr); 2365 { 2366 PetscInt numRootIndices, numRootMatrices; 2367 2368 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 2369 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr); 2370 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr); 2371 for (p = pStartC; p < pEndC; p++) { 2372 PetscInt numRowIndices, numColIndices, matSize, dof; 2373 PetscInt pIndOff, pMatOff, f; 2374 PetscInt *pInd; 2375 PetscInt maxChildId = maxChildIds[p - pStartC]; 2376 PetscScalar *pMat = NULL; 2377 2378 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2379 if (!numColIndices) { 2380 continue; 2381 } 2382 for (f = 0; f <= numFields; f++) { 2383 offsets[f] = 0; 2384 newOffsets[f] = 0; 2385 offsetsCopy[f] = 0; 2386 newOffsetsCopy[f] = 0; 2387 } 2388 numColIndices -= 2 * numFields; 2389 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2390 pInd = &(rootIndices[pIndOff]); 2391 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr); 2392 if (matSize) { 2393 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2394 pMat = &rootMatrices[pMatOff]; 2395 } 2396 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2397 if (dof < 0) { 2398 dof = -(dof + 1); 2399 } 2400 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2401 PetscInt i, j; 2402 PetscInt numRowIndices = matSize / numColIndices; 2403 2404 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2405 PetscInt numIndices, *indices; 2406 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2407 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2408 for (i = 0; i < numColIndices; i++) { 2409 pInd[i] = indices[i]; 2410 } 2411 for (i = 0; i < numFields; i++) { 2412 pInd[numColIndices + i] = offsets[i+1]; 2413 pInd[numColIndices + numFields + i] = offsets[i+1]; 2414 } 2415 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2416 } 2417 else { 2418 PetscInt closureSize, *closure = NULL, cl; 2419 PetscScalar *pMatIn, *pMatModified; 2420 PetscInt numPoints,*points; 2421 2422 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);CHKERRQ(ierr); 2423 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2424 for (j = 0; j < numRowIndices; j++) { 2425 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2426 } 2427 } 2428 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2429 for (f = 0; f < maxFields; f++) { 2430 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2431 else {ierr = PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2432 } 2433 if (numFields) { 2434 for (cl = 0; cl < closureSize; cl++) { 2435 PetscInt c = closure[2 * cl]; 2436 2437 for (f = 0; f < numFields; f++) { 2438 PetscInt fDof; 2439 2440 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr); 2441 offsets[f + 1] += fDof; 2442 } 2443 } 2444 for (f = 0; f < numFields; f++) { 2445 offsets[f + 1] += offsets[f]; 2446 newOffsets[f + 1] = offsets[f + 1]; 2447 } 2448 } 2449 /* TODO : flips here ? */ 2450 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2451 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2452 for (f = 0; f < maxFields; f++) { 2453 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2454 else {ierr = PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2455 } 2456 for (f = 0; f < maxFields; f++) { 2457 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2458 else {ierr = PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2459 } 2460 if (!numFields) { 2461 for (i = 0; i < numRowIndices * numColIndices; i++) { 2462 pMat[i] = pMatModified[i]; 2463 } 2464 } 2465 else { 2466 PetscInt i, j, count; 2467 for (f = 0, count = 0; f < numFields; f++) { 2468 for (i = offsets[f]; i < offsets[f+1]; i++) { 2469 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2470 pMat[count] = pMatModified[i * numColIndices + j]; 2471 } 2472 } 2473 } 2474 } 2475 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified);CHKERRQ(ierr); 2476 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2477 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn);CHKERRQ(ierr); 2478 if (numFields) { 2479 for (f = 0; f < numFields; f++) { 2480 pInd[numColIndices + f] = offsets[f+1]; 2481 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2482 } 2483 for (cl = 0; cl < numPoints; cl++) { 2484 PetscInt globalOff, c = points[2*cl]; 2485 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2486 DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd); 2487 } 2488 } else { 2489 for (cl = 0; cl < numPoints; cl++) { 2490 PetscInt c = points[2*cl], globalOff; 2491 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2492 2493 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2494 DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd); 2495 } 2496 } 2497 for (f = 0; f < maxFields; f++) { 2498 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2499 else {ierr = PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2500 } 2501 ierr = DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);CHKERRQ(ierr); 2502 } 2503 } 2504 else if (matSize) { 2505 PetscInt cOff; 2506 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2507 2508 numRowIndices = matSize / numColIndices; 2509 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2510 ierr = DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2511 ierr = DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);CHKERRQ(ierr); 2512 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr); 2513 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2514 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2515 if (numFields) { 2516 for (f = 0; f < numFields; f++) { 2517 PetscInt fDof; 2518 2519 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr); 2520 offsets[f + 1] = fDof; 2521 for (a = 0; a < aDof; a++) { 2522 PetscInt anchor = anchors[a + aOff]; 2523 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2524 newOffsets[f + 1] += fDof; 2525 } 2526 } 2527 for (f = 0; f < numFields; f++) { 2528 offsets[f + 1] += offsets[f]; 2529 offsetsCopy[f + 1] = offsets[f + 1]; 2530 newOffsets[f + 1] += newOffsets[f]; 2531 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2532 } 2533 ierr = DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);CHKERRQ(ierr); 2534 for (a = 0; a < aDof; a++) { 2535 PetscInt anchor = anchors[a + aOff], lOff; 2536 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2537 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);CHKERRQ(ierr); 2538 } 2539 } 2540 else { 2541 ierr = DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);CHKERRQ(ierr); 2542 for (a = 0; a < aDof; a++) { 2543 PetscInt anchor = anchors[a + aOff], lOff; 2544 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2545 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);CHKERRQ(ierr); 2546 } 2547 } 2548 if (numFields) { 2549 PetscInt count, a; 2550 2551 for (f = 0, count = 0; f < numFields; f++) { 2552 PetscInt iSize = offsets[f + 1] - offsets[f]; 2553 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2554 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2555 count += iSize * jSize; 2556 pInd[numColIndices + f] = offsets[f+1]; 2557 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2558 } 2559 for (a = 0; a < aDof; a++) { 2560 PetscInt anchor = anchors[a + aOff]; 2561 PetscInt gOff; 2562 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2563 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2564 } 2565 } 2566 else { 2567 PetscInt a; 2568 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2569 for (a = 0; a < aDof; a++) { 2570 PetscInt anchor = anchors[a + aOff]; 2571 PetscInt gOff; 2572 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2573 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2574 } 2575 } 2576 ierr = DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);CHKERRQ(ierr); 2577 ierr = DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2578 } 2579 else { 2580 PetscInt gOff; 2581 2582 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 2583 if (numFields) { 2584 for (f = 0; f < numFields; f++) { 2585 PetscInt fDof; 2586 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2587 offsets[f + 1] = fDof + offsets[f]; 2588 } 2589 for (f = 0; f < numFields; f++) { 2590 pInd[numColIndices + f] = offsets[f+1]; 2591 pInd[numColIndices + numFields + f] = offsets[f+1]; 2592 } 2593 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2594 } else { 2595 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2596 } 2597 } 2598 } 2599 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 2600 } 2601 { 2602 PetscSF indicesSF, matricesSF; 2603 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2604 2605 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 2606 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr); 2607 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr); 2608 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr); 2609 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr); 2610 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr); 2611 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 2612 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr); 2613 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr); 2614 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr); 2615 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr); 2616 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr); 2617 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2618 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2619 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2620 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2621 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr); 2622 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 2623 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr); 2624 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 2625 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr); 2626 } 2627 /* count to preallocate */ 2628 ierr = DMGetSection(fine,&localFine);CHKERRQ(ierr); 2629 { 2630 PetscInt nGlobal; 2631 PetscInt *dnnz, *onnz; 2632 PetscLayout rowMap, colMap; 2633 PetscInt rowStart, rowEnd, colStart, colEnd; 2634 PetscInt maxDof; 2635 PetscInt *rowIndices; 2636 DM refTree; 2637 PetscInt **refPointFieldN; 2638 PetscScalar ***refPointFieldMats; 2639 PetscSection refConSec, refAnSec; 2640 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2641 PetscScalar *pointWork; 2642 2643 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr); 2644 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr); 2645 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 2646 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 2647 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 2648 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 2649 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 2650 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 2651 ierr = PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 2652 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2653 for (p = leafStart; p < leafEnd; p++) { 2654 PetscInt gDof, gcDof, gOff; 2655 PetscInt numColIndices, pIndOff, *pInd; 2656 PetscInt matSize; 2657 PetscInt i; 2658 2659 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2660 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2661 if ((gDof - gcDof) <= 0) { 2662 continue; 2663 } 2664 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2665 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2666 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2667 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2668 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2669 numColIndices -= 2 * numFields; 2670 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2671 pInd = &leafIndices[pIndOff]; 2672 offsets[0] = 0; 2673 offsetsCopy[0] = 0; 2674 newOffsets[0] = 0; 2675 newOffsetsCopy[0] = 0; 2676 if (numFields) { 2677 PetscInt f; 2678 for (f = 0; f < numFields; f++) { 2679 PetscInt rowDof; 2680 2681 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2682 offsets[f + 1] = offsets[f] + rowDof; 2683 offsetsCopy[f + 1] = offsets[f + 1]; 2684 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2685 numD[f] = 0; 2686 numO[f] = 0; 2687 } 2688 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2689 for (f = 0; f < numFields; f++) { 2690 PetscInt colOffset = newOffsets[f]; 2691 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2692 2693 for (i = 0; i < numFieldCols; i++) { 2694 PetscInt gInd = pInd[i + colOffset]; 2695 2696 if (gInd >= colStart && gInd < colEnd) { 2697 numD[f]++; 2698 } 2699 else if (gInd >= 0) { /* negative means non-entry */ 2700 numO[f]++; 2701 } 2702 } 2703 } 2704 } 2705 else { 2706 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2707 numD[0] = 0; 2708 numO[0] = 0; 2709 for (i = 0; i < numColIndices; i++) { 2710 PetscInt gInd = pInd[i]; 2711 2712 if (gInd >= colStart && gInd < colEnd) { 2713 numD[0]++; 2714 } 2715 else if (gInd >= 0) { /* negative means non-entry */ 2716 numO[0]++; 2717 } 2718 } 2719 } 2720 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2721 if (!matSize) { /* incoming matrix is identity */ 2722 PetscInt childId; 2723 2724 childId = childIds[p-pStartF]; 2725 if (childId < 0) { /* no child interpolation: one nnz per */ 2726 if (numFields) { 2727 PetscInt f; 2728 for (f = 0; f < numFields; f++) { 2729 PetscInt numRows = offsets[f+1] - offsets[f], row; 2730 for (row = 0; row < numRows; row++) { 2731 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2732 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2733 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2734 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2735 dnnz[gIndFine - rowStart] = 1; 2736 } 2737 else if (gIndCoarse >= 0) { /* remote */ 2738 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2739 onnz[gIndFine - rowStart] = 1; 2740 } 2741 else { /* constrained */ 2742 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2743 } 2744 } 2745 } 2746 } 2747 else { 2748 PetscInt i; 2749 for (i = 0; i < gDof; i++) { 2750 PetscInt gIndCoarse = pInd[i]; 2751 PetscInt gIndFine = rowIndices[i]; 2752 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2753 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2754 dnnz[gIndFine - rowStart] = 1; 2755 } 2756 else if (gIndCoarse >= 0) { /* remote */ 2757 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2758 onnz[gIndFine - rowStart] = 1; 2759 } 2760 else { /* constrained */ 2761 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2762 } 2763 } 2764 } 2765 } 2766 else { /* interpolate from all */ 2767 if (numFields) { 2768 PetscInt f; 2769 for (f = 0; f < numFields; f++) { 2770 PetscInt numRows = offsets[f+1] - offsets[f], row; 2771 for (row = 0; row < numRows; row++) { 2772 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2773 if (gIndFine >= 0) { 2774 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2775 dnnz[gIndFine - rowStart] = numD[f]; 2776 onnz[gIndFine - rowStart] = numO[f]; 2777 } 2778 } 2779 } 2780 } 2781 else { 2782 PetscInt i; 2783 for (i = 0; i < gDof; i++) { 2784 PetscInt gIndFine = rowIndices[i]; 2785 if (gIndFine >= 0) { 2786 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2787 dnnz[gIndFine - rowStart] = numD[0]; 2788 onnz[gIndFine - rowStart] = numO[0]; 2789 } 2790 } 2791 } 2792 } 2793 } 2794 else { /* interpolate from all */ 2795 if (numFields) { 2796 PetscInt f; 2797 for (f = 0; f < numFields; f++) { 2798 PetscInt numRows = offsets[f+1] - offsets[f], row; 2799 for (row = 0; row < numRows; row++) { 2800 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2801 if (gIndFine >= 0) { 2802 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2803 dnnz[gIndFine - rowStart] = numD[f]; 2804 onnz[gIndFine - rowStart] = numO[f]; 2805 } 2806 } 2807 } 2808 } 2809 else { /* every dof get a full row */ 2810 PetscInt i; 2811 for (i = 0; i < gDof; i++) { 2812 PetscInt gIndFine = rowIndices[i]; 2813 if (gIndFine >= 0) { 2814 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2815 dnnz[gIndFine - rowStart] = numD[0]; 2816 onnz[gIndFine - rowStart] = numO[0]; 2817 } 2818 } 2819 } 2820 } 2821 } 2822 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr); 2823 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2824 2825 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 2826 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2827 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 2828 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 2829 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 2830 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr); 2831 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr); 2832 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr); 2833 for (p = leafStart; p < leafEnd; p++) { 2834 PetscInt gDof, gcDof, gOff; 2835 PetscInt numColIndices, pIndOff, *pInd; 2836 PetscInt matSize; 2837 PetscInt childId; 2838 2839 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2840 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2841 if ((gDof - gcDof) <= 0) { 2842 continue; 2843 } 2844 childId = childIds[p-pStartF]; 2845 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2846 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2847 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2848 numColIndices -= 2 * numFields; 2849 pInd = &leafIndices[pIndOff]; 2850 offsets[0] = 0; 2851 offsetsCopy[0] = 0; 2852 newOffsets[0] = 0; 2853 newOffsetsCopy[0] = 0; 2854 rowOffsets[0] = 0; 2855 if (numFields) { 2856 PetscInt f; 2857 for (f = 0; f < numFields; f++) { 2858 PetscInt rowDof; 2859 2860 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2861 offsets[f + 1] = offsets[f] + rowDof; 2862 offsetsCopy[f + 1] = offsets[f + 1]; 2863 rowOffsets[f + 1] = pInd[numColIndices + f]; 2864 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2865 } 2866 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2867 } 2868 else { 2869 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2870 } 2871 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2872 if (!matSize) { /* incoming matrix is identity */ 2873 if (childId < 0) { /* no child interpolation: scatter */ 2874 if (numFields) { 2875 PetscInt f; 2876 for (f = 0; f < numFields; f++) { 2877 PetscInt numRows = offsets[f+1] - offsets[f], row; 2878 for (row = 0; row < numRows; row++) { 2879 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2880 } 2881 } 2882 } 2883 else { 2884 PetscInt numRows = gDof, row; 2885 for (row = 0; row < numRows; row++) { 2886 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2887 } 2888 } 2889 } 2890 else { /* interpolate from all */ 2891 if (numFields) { 2892 PetscInt f; 2893 for (f = 0; f < numFields; f++) { 2894 PetscInt numRows = offsets[f+1] - offsets[f]; 2895 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2896 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2897 } 2898 } 2899 else { 2900 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2901 } 2902 } 2903 } 2904 else { /* interpolate from all */ 2905 PetscInt pMatOff; 2906 PetscScalar *pMat; 2907 2908 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2909 pMat = &leafMatrices[pMatOff]; 2910 if (childId < 0) { /* copy the incoming matrix */ 2911 if (numFields) { 2912 PetscInt f, count; 2913 for (f = 0, count = 0; f < numFields; f++) { 2914 PetscInt numRows = offsets[f+1]-offsets[f]; 2915 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2916 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2917 PetscScalar *inMat = &pMat[count]; 2918 2919 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2920 count += numCols * numInRows; 2921 } 2922 } 2923 else { 2924 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2925 } 2926 } 2927 else { /* multiply the incoming matrix by the child interpolation */ 2928 if (numFields) { 2929 PetscInt f, count; 2930 for (f = 0, count = 0; f < numFields; f++) { 2931 PetscInt numRows = offsets[f+1]-offsets[f]; 2932 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2933 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2934 PetscScalar *inMat = &pMat[count]; 2935 PetscInt i, j, k; 2936 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2937 for (i = 0; i < numRows; i++) { 2938 for (j = 0; j < numCols; j++) { 2939 PetscScalar val = 0.; 2940 for (k = 0; k < numInRows; k++) { 2941 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2942 } 2943 pointWork[i * numCols + j] = val; 2944 } 2945 } 2946 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2947 count += numCols * numInRows; 2948 } 2949 } 2950 else { /* every dof gets a full row */ 2951 PetscInt numRows = gDof; 2952 PetscInt numCols = numColIndices; 2953 PetscInt numInRows = matSize / numColIndices; 2954 PetscInt i, j, k; 2955 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2956 for (i = 0; i < numRows; i++) { 2957 for (j = 0; j < numCols; j++) { 2958 PetscScalar val = 0.; 2959 for (k = 0; k < numInRows; k++) { 2960 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2961 } 2962 pointWork[i * numCols + j] = val; 2963 } 2964 } 2965 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 2966 } 2967 } 2968 } 2969 } 2970 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2971 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2972 ierr = PetscFree(pointWork);CHKERRQ(ierr); 2973 } 2974 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2975 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2976 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 2977 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 2978 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 2979 ierr = PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);CHKERRQ(ierr); 2980 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 2981 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 2982 PetscFunctionReturn(0); 2983 } 2984 2985 /* 2986 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2987 * 2988 * for each coarse dof \phi^c_i: 2989 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2990 * for each fine dof \phi^f_j; 2991 * a_{i,j} = 0; 2992 * for each fine dof \phi^f_k: 2993 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2994 * [^^^ this is = \phi^c_i ^^^] 2995 */ 2996 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2997 { 2998 PetscDS ds; 2999 PetscSection section, cSection; 3000 DMLabel canonical, depth; 3001 Mat cMat, mat; 3002 PetscInt *nnz; 3003 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 3004 PetscInt m, n; 3005 PetscScalar *pointScalar; 3006 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 3007 PetscErrorCode ierr; 3008 3009 PetscFunctionBegin; 3010 ierr = DMGetSection(refTree,§ion);CHKERRQ(ierr); 3011 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr); 3012 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr); 3013 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr); 3014 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3015 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3016 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr); 3017 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr); 3018 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr); 3019 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr); 3020 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr); 3021 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr); 3022 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */ 3023 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3024 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 3025 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3026 const PetscInt *children; 3027 PetscInt numChildren; 3028 PetscInt i, numChildDof, numSelfDof; 3029 3030 if (canonical) { 3031 PetscInt pCanonical; 3032 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3033 if (p != pCanonical) continue; 3034 } 3035 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3036 if (!numChildren) continue; 3037 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3038 PetscInt child = children[i]; 3039 PetscInt dof; 3040 3041 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3042 numChildDof += dof; 3043 } 3044 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3045 if (!numChildDof || !numSelfDof) continue; 3046 for (f = 0; f < numFields; f++) { 3047 PetscInt selfOff; 3048 3049 if (numSecFields) { /* count the dofs for just this field */ 3050 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3051 PetscInt child = children[i]; 3052 PetscInt dof; 3053 3054 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3055 numChildDof += dof; 3056 } 3057 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3058 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3059 } 3060 else { 3061 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3062 } 3063 for (i = 0; i < numSelfDof; i++) { 3064 nnz[selfOff + i] = numChildDof; 3065 } 3066 } 3067 } 3068 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr); 3069 ierr = PetscFree(nnz);CHKERRQ(ierr); 3070 /* Setp 2: compute entries */ 3071 for (p = pStart; p < pEnd; p++) { 3072 const PetscInt *children; 3073 PetscInt numChildren; 3074 PetscInt i, numChildDof, numSelfDof; 3075 3076 /* same conditions about when entries occur */ 3077 if (canonical) { 3078 PetscInt pCanonical; 3079 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3080 if (p != pCanonical) continue; 3081 } 3082 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3083 if (!numChildren) continue; 3084 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3085 PetscInt child = children[i]; 3086 PetscInt dof; 3087 3088 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3089 numChildDof += dof; 3090 } 3091 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3092 if (!numChildDof || !numSelfDof) continue; 3093 3094 for (f = 0; f < numFields; f++) { 3095 PetscInt selfOff, Nc, parentCell; 3096 PetscInt cellShapeOff; 3097 PetscObject disc; 3098 PetscDualSpace dsp; 3099 PetscClassId classId; 3100 PetscScalar *pointMat; 3101 PetscInt *matRows, *matCols; 3102 PetscInt pO = PETSC_MIN_INT; 3103 const PetscInt *depthNumDof; 3104 3105 if (numSecFields) { 3106 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3107 PetscInt child = children[i]; 3108 PetscInt dof; 3109 3110 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3111 numChildDof += dof; 3112 } 3113 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3114 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3115 } 3116 else { 3117 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3118 } 3119 3120 /* find a cell whose closure contains p */ 3121 if (p >= cStart && p < cEnd) { 3122 parentCell = p; 3123 } 3124 else { 3125 PetscInt *star = NULL; 3126 PetscInt numStar; 3127 3128 parentCell = -1; 3129 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3130 for (i = numStar - 1; i >= 0; i--) { 3131 PetscInt c = star[2 * i]; 3132 3133 if (c >= cStart && c < cEnd) { 3134 parentCell = c; 3135 break; 3136 } 3137 } 3138 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3139 } 3140 /* determine the offset of p's shape functions withing parentCell's shape functions */ 3141 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3142 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr); 3143 if (classId == PETSCFE_CLASSID) { 3144 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr); 3145 } 3146 else if (classId == PETSCFV_CLASSID) { 3147 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr); 3148 } 3149 else { 3150 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object"); 3151 } 3152 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr); 3153 ierr = PetscDualSpaceGetNumComponents(dsp,&Nc);CHKERRQ(ierr); 3154 { 3155 PetscInt *closure = NULL; 3156 PetscInt numClosure; 3157 3158 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3159 for (i = 0, cellShapeOff = 0; i < numClosure; i++) { 3160 PetscInt point = closure[2 * i], pointDepth; 3161 3162 pO = closure[2 * i + 1]; 3163 if (point == p) break; 3164 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3165 cellShapeOff += depthNumDof[pointDepth]; 3166 } 3167 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3168 } 3169 3170 ierr = DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr); 3171 ierr = DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr); 3172 matCols = matRows + numSelfDof; 3173 for (i = 0; i < numSelfDof; i++) { 3174 matRows[i] = selfOff + i; 3175 } 3176 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.; 3177 { 3178 PetscInt colOff = 0; 3179 3180 for (i = 0; i < numChildren; i++) { 3181 PetscInt child = children[i]; 3182 PetscInt dof, off, j; 3183 3184 if (numSecFields) { 3185 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr); 3186 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr); 3187 } 3188 else { 3189 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr); 3190 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr); 3191 } 3192 3193 for (j = 0; j < dof; j++) { 3194 matCols[colOff++] = off + j; 3195 } 3196 } 3197 } 3198 if (classId == PETSCFE_CLASSID) { 3199 PetscFE fe = (PetscFE) disc; 3200 PetscInt fSize; 3201 3202 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 3203 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr); 3204 for (i = 0; i < numSelfDof; i++) { /* for every shape function */ 3205 PetscQuadrature q; 3206 PetscInt dim, thisNc, numPoints, j, k; 3207 const PetscReal *points; 3208 const PetscReal *weights; 3209 PetscInt *closure = NULL; 3210 PetscInt numClosure; 3211 PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfDof - 1 - i) : i); 3212 PetscReal *Bparent; 3213 3214 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr); 3215 ierr = PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);CHKERRQ(ierr); 3216 if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc); 3217 ierr = PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3218 for (j = 0; j < numPoints; j++) { 3219 PetscInt childCell = -1; 3220 PetscReal *parentValAtPoint; 3221 const PetscReal xi0[3] = {-1.,-1.,-1.}; 3222 const PetscReal *pointReal = &points[dim * j]; 3223 const PetscScalar *point; 3224 PetscReal *Bchild; 3225 PetscInt childCellShapeOff, pointMatOff; 3226 #if defined(PETSC_USE_COMPLEX) 3227 PetscInt d; 3228 3229 for (d = 0; d < dim; d++) { 3230 pointScalar[d] = points[dim * j + d]; 3231 } 3232 point = pointScalar; 3233 #else 3234 point = pointReal; 3235 #endif 3236 3237 parentValAtPoint = &Bparent[(fSize * j + parentCellShapeDof) * Nc]; 3238 3239 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3240 PetscInt child = children[k]; 3241 PetscInt *star = NULL; 3242 PetscInt numStar, s; 3243 3244 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3245 for (s = numStar - 1; s >= 0; s--) { 3246 PetscInt c = star[2 * s]; 3247 3248 if (c < cStart || c >= cEnd) continue; 3249 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr); 3250 if (childCell >= 0) break; 3251 } 3252 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3253 if (childCell >= 0) break; 3254 } 3255 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point"); 3256 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 3257 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr); 3258 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp); 3259 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef); 3260 3261 ierr = PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3262 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3263 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3264 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3265 PetscInt l; 3266 3267 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr); 3268 childDof = depthNumDof[childDepth]; 3269 for (l = 0, childCellShapeOff = 0; l < numClosure; l++) { 3270 PetscInt point = closure[2 * l]; 3271 PetscInt pointDepth; 3272 3273 childO = closure[2 * l + 1]; 3274 if (point == child) break; 3275 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3276 childCellShapeOff += depthNumDof[pointDepth]; 3277 } 3278 if (l == numClosure) { 3279 pointMatOff += childDof; 3280 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3281 } 3282 for (l = 0; l < childDof; l++) { 3283 PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l), m; 3284 PetscReal *childValAtPoint; 3285 PetscReal val = 0.; 3286 3287 childValAtPoint = &Bchild[childCellDof * Nc]; 3288 for (m = 0; m < Nc; m++) { 3289 val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m]; 3290 } 3291 3292 pointMat[i * numChildDof + pointMatOff + l] += val; 3293 } 3294 pointMatOff += childDof; 3295 } 3296 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3297 ierr = PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3298 } 3299 ierr = PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); 3300 } 3301 } 3302 else { /* just the volume-weighted averages of the children */ 3303 PetscReal parentVol; 3304 PetscInt childCell; 3305 3306 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr); 3307 for (i = 0, childCell = 0; i < numChildren; i++) { 3308 PetscInt child = children[i], j; 3309 PetscReal childVol; 3310 3311 if (child < cStart || child >= cEnd) continue; 3312 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr); 3313 for (j = 0; j < Nc; j++) { 3314 pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol; 3315 } 3316 childCell++; 3317 } 3318 } 3319 /* Insert pointMat into mat */ 3320 ierr = MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr); 3321 ierr = DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr); 3322 ierr = DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr); 3323 } 3324 } 3325 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr); 3326 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr); 3327 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3328 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3329 *inj = mat; 3330 PetscFunctionReturn(0); 3331 } 3332 3333 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3334 { 3335 PetscDS ds; 3336 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3337 PetscScalar ***refPointFieldMats; 3338 PetscSection refConSec, refSection; 3339 PetscErrorCode ierr; 3340 3341 PetscFunctionBegin; 3342 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3343 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3344 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3345 ierr = DMGetSection(refTree,&refSection);CHKERRQ(ierr); 3346 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3347 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 3348 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 3349 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 3350 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr); 3351 for (p = pRefStart; p < pRefEnd; p++) { 3352 PetscInt parent, pDof, parentDof; 3353 3354 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3355 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3356 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3357 if (!pDof || !parentDof || parent == p) continue; 3358 3359 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 3360 for (f = 0; f < numFields; f++) { 3361 PetscInt cDof, cOff, numCols, r; 3362 3363 if (numFields > 1) { 3364 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3365 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 3366 } 3367 else { 3368 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3369 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 3370 } 3371 3372 for (r = 0; r < cDof; r++) { 3373 rows[r] = cOff + r; 3374 } 3375 numCols = 0; 3376 { 3377 PetscInt aDof, aOff, j; 3378 3379 if (numFields > 1) { 3380 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr); 3381 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr); 3382 } 3383 else { 3384 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr); 3385 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr); 3386 } 3387 3388 for (j = 0; j < aDof; j++) { 3389 cols[numCols++] = aOff + j; 3390 } 3391 } 3392 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3393 /* transpose of constraint matrix */ 3394 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3395 } 3396 } 3397 *childrenMats = refPointFieldMats; 3398 ierr = PetscFree(rows);CHKERRQ(ierr); 3399 ierr = PetscFree(cols);CHKERRQ(ierr); 3400 PetscFunctionReturn(0); 3401 } 3402 3403 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3404 { 3405 PetscDS ds; 3406 PetscScalar ***refPointFieldMats; 3407 PetscInt numFields, pRefStart, pRefEnd, p, f; 3408 PetscSection refConSec, refSection; 3409 PetscErrorCode ierr; 3410 3411 PetscFunctionBegin; 3412 refPointFieldMats = *childrenMats; 3413 *childrenMats = NULL; 3414 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3415 ierr = DMGetSection(refTree,&refSection);CHKERRQ(ierr); 3416 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3417 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3418 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3419 for (p = pRefStart; p < pRefEnd; p++) { 3420 PetscInt parent, pDof, parentDof; 3421 3422 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3423 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3424 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3425 if (!pDof || !parentDof || parent == p) continue; 3426 3427 for (f = 0; f < numFields; f++) { 3428 PetscInt cDof; 3429 3430 if (numFields > 1) { 3431 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3432 } 3433 else { 3434 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3435 } 3436 3437 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 3438 } 3439 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 3440 } 3441 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 3442 PetscFunctionReturn(0); 3443 } 3444 3445 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3446 { 3447 Mat cMatRef; 3448 PetscObject injRefObj; 3449 PetscErrorCode ierr; 3450 3451 PetscFunctionBegin; 3452 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr); 3453 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr); 3454 *injRef = (Mat) injRefObj; 3455 if (!*injRef) { 3456 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr); 3457 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr); 3458 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3459 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr); 3460 } 3461 PetscFunctionReturn(0); 3462 } 3463 3464 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) 3465 { 3466 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3467 PetscSection globalCoarse, globalFine; 3468 PetscSection localCoarse, localFine, leafIndicesSec; 3469 PetscSection multiRootSec, rootIndicesSec; 3470 PetscInt *leafInds, *rootInds = NULL; 3471 const PetscInt *rootDegrees; 3472 PetscScalar *leafVals = NULL, *rootVals = NULL; 3473 PetscSF coarseToFineEmbedded; 3474 PetscErrorCode ierr; 3475 3476 PetscFunctionBegin; 3477 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3478 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3479 ierr = DMGetSection(fine,&localFine);CHKERRQ(ierr); 3480 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3481 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 3482 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr); 3483 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr); 3484 { /* winnow fine points that don't have global dofs out of the sf */ 3485 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3486 const PetscInt *leaves; 3487 3488 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3489 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3490 p = leaves ? leaves[l] : l; 3491 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3492 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3493 if ((dof - cdof) > 0) { 3494 numPointsWithDofs++; 3495 3496 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr); 3497 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr); 3498 } 3499 } 3500 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3501 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr); 3502 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr); 3503 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr); 3504 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);} 3505 for (l = 0, offset = 0; l < nleaves; l++) { 3506 p = leaves ? leaves[l] : l; 3507 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3508 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3509 if ((dof - cdof) > 0) { 3510 PetscInt off, gOff; 3511 PetscInt *pInd; 3512 PetscScalar *pVal = NULL; 3513 3514 pointsWithDofs[offset++] = l; 3515 3516 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3517 3518 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3519 if (gatheredValues) { 3520 PetscInt i; 3521 3522 pVal = &leafVals[off + 1]; 3523 for (i = 0; i < dof; i++) pVal[i] = 0.; 3524 } 3525 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 3526 3527 offsets[0] = 0; 3528 if (numFields) { 3529 PetscInt f; 3530 3531 for (f = 0; f < numFields; f++) { 3532 PetscInt fDof; 3533 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr); 3534 offsets[f + 1] = fDof + offsets[f]; 3535 } 3536 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 3537 } else { 3538 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 3539 } 3540 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);} 3541 } 3542 } 3543 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 3544 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3545 } 3546 3547 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3548 ierr = DMGetSection(coarse,&localCoarse);CHKERRQ(ierr); 3549 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3550 3551 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3552 MPI_Datatype threeInt; 3553 PetscMPIInt rank; 3554 PetscInt (*parentNodeAndIdCoarse)[3]; 3555 PetscInt (*parentNodeAndIdFine)[3]; 3556 PetscInt p, nleaves, nleavesToParents; 3557 PetscSF pointSF, sfToParents; 3558 const PetscInt *ilocal; 3559 const PetscSFNode *iremote; 3560 PetscSFNode *iremoteToParents; 3561 PetscInt *ilocalToParents; 3562 3563 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr); 3564 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr); 3565 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr); 3566 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr); 3567 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr); 3568 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 3569 for (p = pStartC; p < pEndC; p++) { 3570 PetscInt parent, childId; 3571 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr); 3572 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3573 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3574 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3575 if (nleaves > 0) { 3576 PetscInt leaf = -1; 3577 3578 if (ilocal) { 3579 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr); 3580 } 3581 else { 3582 leaf = p - pStartC; 3583 } 3584 if (leaf >= 0) { 3585 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3586 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3587 } 3588 } 3589 } 3590 for (p = pStartF; p < pEndF; p++) { 3591 parentNodeAndIdFine[p - pStartF][0] = -1; 3592 parentNodeAndIdFine[p - pStartF][1] = -1; 3593 parentNodeAndIdFine[p - pStartF][2] = -1; 3594 } 3595 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3596 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3597 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3598 PetscInt dof; 3599 3600 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr); 3601 if (dof) { 3602 PetscInt off; 3603 3604 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3605 if (gatheredIndices) { 3606 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3607 } else if (gatheredValues) { 3608 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3609 } 3610 } 3611 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3612 nleavesToParents++; 3613 } 3614 } 3615 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr); 3616 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr); 3617 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3618 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3619 ilocalToParents[nleavesToParents] = p - pStartF; 3620 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3621 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3622 nleavesToParents++; 3623 } 3624 } 3625 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr); 3626 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr); 3627 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3628 3629 coarseToFineEmbedded = sfToParents; 3630 3631 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3632 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr); 3633 } 3634 3635 { /* winnow out coarse points that don't have dofs */ 3636 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3637 PetscSF sfDofsOnly; 3638 3639 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3640 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3641 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3642 if ((dof - cdof) > 0) { 3643 numPointsWithDofs++; 3644 } 3645 } 3646 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3647 for (p = pStartC, offset = 0; p < pEndC; p++) { 3648 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3649 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3650 if ((dof - cdof) > 0) { 3651 pointsWithDofs[offset++] = p - pStartC; 3652 } 3653 } 3654 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr); 3655 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3656 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3657 coarseToFineEmbedded = sfDofsOnly; 3658 } 3659 3660 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3661 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3662 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3663 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr); 3664 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr); 3665 for (p = pStartC; p < pEndC; p++) { 3666 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr); 3667 } 3668 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr); 3669 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr); 3670 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 3671 { /* distribute the leaf section */ 3672 PetscSF multi, multiInv, indicesSF; 3673 PetscInt *remoteOffsets, numRootIndices; 3674 3675 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr); 3676 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr); 3677 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr); 3678 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr); 3679 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 3680 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr); 3681 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 3682 if (gatheredIndices) { 3683 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr); 3684 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3685 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3686 } 3687 if (gatheredValues) { 3688 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr); 3689 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3690 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3691 } 3692 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 3693 } 3694 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3695 ierr = PetscFree(leafInds);CHKERRQ(ierr); 3696 ierr = PetscFree(leafVals);CHKERRQ(ierr); 3697 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3698 *rootMultiSec = multiRootSec; 3699 *multiLeafSec = rootIndicesSec; 3700 if (gatheredIndices) *gatheredIndices = rootInds; 3701 if (gatheredValues) *gatheredValues = rootVals; 3702 PetscFunctionReturn(0); 3703 } 3704 3705 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3706 { 3707 DM refTree; 3708 PetscSection multiRootSec, rootIndicesSec; 3709 PetscSection globalCoarse, globalFine; 3710 PetscSection localCoarse, localFine; 3711 PetscSection cSecRef; 3712 PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd; 3713 Mat injRef; 3714 PetscInt numFields, maxDof; 3715 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3716 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3717 PetscLayout rowMap, colMap; 3718 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3719 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3720 PetscErrorCode ierr; 3721 3722 PetscFunctionBegin; 3723 3724 /* get the templates for the fine-to-coarse injection from the reference tree */ 3725 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 3726 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 3727 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3728 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 3729 3730 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3731 ierr = DMGetSection(fine,&localFine);CHKERRQ(ierr); 3732 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3733 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 3734 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3735 ierr = DMGetSection(coarse,&localCoarse);CHKERRQ(ierr); 3736 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3737 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 3738 { 3739 PetscInt maxFields = PetscMax(1,numFields) + 1; 3740 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 3741 } 3742 3743 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr); 3744 3745 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr); 3746 3747 /* count indices */ 3748 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 3749 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 3750 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 3751 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 3752 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 3753 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr); 3754 for (p = pStartC; p < pEndC; p++) { 3755 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3756 3757 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3758 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3759 if ((dof - cdof) <= 0) continue; 3760 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3761 3762 rowOffsets[0] = 0; 3763 offsetsCopy[0] = 0; 3764 if (numFields) { 3765 PetscInt f; 3766 3767 for (f = 0; f < numFields; f++) { 3768 PetscInt fDof; 3769 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3770 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3771 } 3772 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3773 } else { 3774 ierr = 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 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3881 } else { 3882 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3883 rowOffsets[1] = offsetsCopy[0]; 3884 } 3885 3886 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3887 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3888 leafEnd = leafStart + numLeaves; 3889 for (l = leafStart; l < leafEnd; l++) { 3890 PetscInt numIndices, childId, offset; 3891 const PetscInt *childIndices; 3892 3893 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3894 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3895 childId = rootIndices[offset++]; 3896 childIndices = &rootIndices[offset]; 3897 numIndices--; 3898 3899 if (childId == -1) { /* equivalent points: scatter */ 3900 PetscInt i; 3901 3902 for (i = 0; i < numIndices; i++) { 3903 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr); 3904 } 3905 } 3906 else { 3907 PetscInt parentId, f, lim; 3908 3909 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3910 3911 lim = PetscMax(1,numFields); 3912 offsets[0] = 0; 3913 if (numFields) { 3914 PetscInt f; 3915 3916 for (f = 0; f < numFields; f++) { 3917 PetscInt fDof; 3918 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3919 3920 offsets[f + 1] = fDof + offsets[f]; 3921 } 3922 } 3923 else { 3924 PetscInt cDof; 3925 3926 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3927 offsets[1] = cDof; 3928 } 3929 for (f = 0; f < lim; f++) { 3930 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3931 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3932 const PetscInt *colIndices = &childIndices[offsets[f]]; 3933 3934 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr); 3935 } 3936 } 3937 } 3938 } 3939 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 3940 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 3941 ierr = PetscFree(parentIndices);CHKERRQ(ierr); 3942 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3943 ierr = PetscFree(rootIndices);CHKERRQ(ierr); 3944 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 3945 3946 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3947 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3948 PetscFunctionReturn(0); 3949 } 3950 3951 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3952 { 3953 PetscDS ds; 3954 PetscSF coarseToFineEmbedded; 3955 PetscSection globalCoarse, globalFine; 3956 PetscSection localCoarse, localFine; 3957 PetscSection aSec, cSec; 3958 PetscSection rootValuesSec; 3959 PetscSection leafValuesSec; 3960 PetscScalar *rootValues, *leafValues; 3961 IS aIS; 3962 const PetscInt *anchors; 3963 Mat cMat; 3964 PetscInt numFields; 3965 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior; 3966 PetscInt aStart, aEnd, cStart, cEnd; 3967 PetscInt *maxChildIds; 3968 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3969 PetscFV fv = NULL; 3970 PetscInt dim, numFVcomps = -1, fvField = -1; 3971 DM cellDM = NULL, gradDM = NULL; 3972 const PetscScalar *cellGeomArray = NULL; 3973 const PetscScalar *gradArray = NULL; 3974 PetscErrorCode ierr; 3975 3976 PetscFunctionBegin; 3977 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 3978 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3979 ierr = DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr); 3980 ierr = DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 3981 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 3982 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3983 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3984 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr); 3985 { /* winnow fine points that don't have global dofs out of the sf */ 3986 PetscInt nleaves, l; 3987 const PetscInt *leaves; 3988 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3989 3990 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3991 3992 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3993 PetscInt p = leaves ? leaves[l] : l; 3994 3995 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3996 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3997 if ((dof - cdof) > 0) { 3998 numPointsWithDofs++; 3999 } 4000 } 4001 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 4002 for (l = 0, offset = 0; l < nleaves; l++) { 4003 PetscInt p = leaves ? leaves[l] : l; 4004 4005 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4006 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4007 if ((dof - cdof) > 0) { 4008 pointsWithDofs[offset++] = l; 4009 } 4010 } 4011 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 4012 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 4013 } 4014 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4015 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 4016 for (p = pStartC; p < pEndC; p++) { 4017 maxChildIds[p - pStartC] = -2; 4018 } 4019 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4020 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4021 4022 ierr = DMGetSection(coarse,&localCoarse);CHKERRQ(ierr); 4023 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4024 4025 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 4026 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 4027 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 4028 4029 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 4030 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 4031 4032 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4033 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr); 4034 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr); 4035 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 4036 { 4037 PetscInt maxFields = PetscMax(1,numFields) + 1; 4038 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 4039 } 4040 if (grad) { 4041 PetscInt i; 4042 4043 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr); 4044 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr); 4045 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr); 4046 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr); 4047 for (i = 0; i < PetscMax(1,numFields); i++) { 4048 PetscObject obj; 4049 PetscClassId id; 4050 4051 ierr = DMGetField(coarse, i, &obj);CHKERRQ(ierr); 4052 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr); 4053 if (id == PETSCFV_CLASSID) { 4054 fv = (PetscFV) obj; 4055 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr); 4056 fvField = i; 4057 break; 4058 } 4059 } 4060 } 4061 4062 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4063 PetscInt dof; 4064 PetscInt maxChildId = maxChildIds[p - pStartC]; 4065 PetscInt numValues = 0; 4066 4067 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4068 if (dof < 0) { 4069 dof = -(dof + 1); 4070 } 4071 offsets[0] = 0; 4072 newOffsets[0] = 0; 4073 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4074 PetscInt *closure = NULL, closureSize, cl; 4075 4076 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4077 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4078 PetscInt c = closure[2 * cl], clDof; 4079 4080 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 4081 numValues += clDof; 4082 } 4083 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4084 } 4085 else if (maxChildId == -1) { 4086 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr); 4087 } 4088 /* we will pack the column indices with the field offsets */ 4089 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4090 /* also send the centroid, and the gradient */ 4091 numValues += dim * (1 + numFVcomps); 4092 } 4093 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr); 4094 } 4095 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr); 4096 { 4097 PetscInt numRootValues; 4098 const PetscScalar *coarseArray; 4099 4100 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr); 4101 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr); 4102 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4103 for (p = pStartC; p < pEndC; p++) { 4104 PetscInt numValues; 4105 PetscInt pValOff; 4106 PetscScalar *pVal; 4107 PetscInt maxChildId = maxChildIds[p - pStartC]; 4108 4109 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr); 4110 if (!numValues) { 4111 continue; 4112 } 4113 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr); 4114 pVal = &(rootValues[pValOff]); 4115 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4116 PetscInt closureSize = numValues; 4117 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr); 4118 if (grad && p >= cellStart && p < cellEnd) { 4119 PetscFVCellGeom *cg; 4120 PetscScalar *gradVals = NULL; 4121 PetscInt i; 4122 4123 pVal += (numValues - dim * (1 + numFVcomps)); 4124 4125 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr); 4126 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4127 pVal += dim; 4128 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr); 4129 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4130 } 4131 } 4132 else if (maxChildId == -1) { 4133 PetscInt lDof, lOff, i; 4134 4135 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr); 4136 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr); 4137 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4138 } 4139 } 4140 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4141 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 4142 } 4143 { 4144 PetscSF valuesSF; 4145 PetscInt *remoteOffsetsValues, numLeafValues; 4146 4147 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr); 4148 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr); 4149 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr); 4150 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 4151 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr); 4152 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr); 4153 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr); 4154 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4155 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4156 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr); 4157 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4158 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr); 4159 } 4160 ierr = DMGetSection(fine,&localFine);CHKERRQ(ierr); 4161 { 4162 PetscInt maxDof; 4163 PetscInt *rowIndices; 4164 DM refTree; 4165 PetscInt **refPointFieldN; 4166 PetscScalar ***refPointFieldMats; 4167 PetscSection refConSec, refAnSec; 4168 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4169 PetscScalar *pointWork; 4170 4171 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 4172 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 4173 ierr = DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr); 4174 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 4175 ierr = DMGetDS(fine,&ds);CHKERRQ(ierr); 4176 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4177 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4178 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 4179 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 4180 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4181 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 4182 ierr = DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4183 ierr = DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4184 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4185 for (p = leafStart; p < leafEnd; p++) { 4186 PetscInt gDof, gcDof, gOff, lDof; 4187 PetscInt numValues, pValOff; 4188 PetscInt childId; 4189 const PetscScalar *pVal; 4190 const PetscScalar *fvGradData = NULL; 4191 4192 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 4193 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr); 4194 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 4195 if ((gDof - gcDof) <= 0) { 4196 continue; 4197 } 4198 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 4199 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr); 4200 if (!numValues) continue; 4201 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr); 4202 pVal = &leafValues[pValOff]; 4203 offsets[0] = 0; 4204 offsetsCopy[0] = 0; 4205 newOffsets[0] = 0; 4206 newOffsetsCopy[0] = 0; 4207 childId = cids[p - pStartF]; 4208 if (numFields) { 4209 PetscInt f; 4210 for (f = 0; f < numFields; f++) { 4211 PetscInt rowDof; 4212 4213 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 4214 offsets[f + 1] = offsets[f] + rowDof; 4215 offsetsCopy[f + 1] = offsets[f + 1]; 4216 /* TODO: closure indices */ 4217 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4218 } 4219 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 4220 } 4221 else { 4222 offsets[0] = 0; 4223 offsets[1] = lDof; 4224 newOffsets[0] = 0; 4225 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4226 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 4227 } 4228 if (childId == -1) { /* no child interpolation: one nnz per */ 4229 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr); 4230 } else { 4231 PetscInt f; 4232 4233 if (grad && p >= cellStart && p < cellEnd) { 4234 numValues -= (dim * (1 + numFVcomps)); 4235 fvGradData = &pVal[numValues]; 4236 } 4237 for (f = 0; f < PetscMax(1,numFields); f++) { 4238 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4239 PetscInt numRows = offsets[f+1] - offsets[f]; 4240 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4241 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4242 PetscScalar *rVal = &pointWork[offsets[f]]; 4243 PetscInt i, j; 4244 4245 #if 0 4246 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); 4247 #endif 4248 for (i = 0; i < numRows; i++) { 4249 PetscScalar val = 0.; 4250 for (j = 0; j < numCols; j++) { 4251 val += childMat[i * numCols + j] * cVal[j]; 4252 } 4253 rVal[i] = val; 4254 } 4255 if (f == fvField && p >= cellStart && p < cellEnd) { 4256 PetscReal centroid[3]; 4257 PetscScalar diff[3]; 4258 const PetscScalar *parentCentroid = &fvGradData[0]; 4259 const PetscScalar *gradient = &fvGradData[dim]; 4260 4261 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr); 4262 for (i = 0; i < dim; i++) { 4263 diff[i] = centroid[i] - parentCentroid[i]; 4264 } 4265 for (i = 0; i < numFVcomps; i++) { 4266 PetscScalar val = 0.; 4267 4268 for (j = 0; j < dim; j++) { 4269 val += gradient[dim * i + j] * diff[j]; 4270 } 4271 rVal[i] += val; 4272 } 4273 } 4274 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr); 4275 } 4276 } 4277 } 4278 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4279 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 4280 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr); 4281 } 4282 ierr = PetscFree(leafValues);CHKERRQ(ierr); 4283 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr); 4284 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 4285 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 4286 PetscFunctionReturn(0); 4287 } 4288 4289 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4290 { 4291 PetscDS ds; 4292 DM refTree; 4293 PetscSection multiRootSec, rootIndicesSec; 4294 PetscSection globalCoarse, globalFine; 4295 PetscSection localCoarse, localFine; 4296 PetscSection cSecRef; 4297 PetscInt *parentIndices, pRefStart, pRefEnd; 4298 PetscScalar *rootValues, *parentValues; 4299 Mat injRef; 4300 PetscInt numFields, maxDof; 4301 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4302 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4303 PetscLayout rowMap, colMap; 4304 PetscInt rowStart, rowEnd, colStart, colEnd; 4305 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4306 PetscErrorCode ierr; 4307 4308 PetscFunctionBegin; 4309 4310 /* get the templates for the fine-to-coarse injection from the reference tree */ 4311 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4312 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4313 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 4314 ierr = DMGetDS(coarse,&ds);CHKERRQ(ierr); 4315 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4316 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 4317 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4318 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 4319 4320 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4321 ierr = DMGetSection(fine,&localFine);CHKERRQ(ierr); 4322 ierr = DMGetGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4323 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 4324 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4325 ierr = DMGetSection(coarse,&localCoarse);CHKERRQ(ierr); 4326 ierr = DMGetGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4327 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 4328 { 4329 PetscInt maxFields = PetscMax(1,numFields) + 1; 4330 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 4331 } 4332 4333 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr); 4334 4335 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr); 4336 4337 /* count indices */ 4338 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr); 4339 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr); 4340 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 4341 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 4342 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 4343 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 4344 /* insert values */ 4345 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4346 for (p = pStartC; p < pEndC; p++) { 4347 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4348 PetscBool contribute = PETSC_FALSE; 4349 4350 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4351 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 4352 if ((dof - cdof) <= 0) continue; 4353 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr); 4354 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 4355 4356 rowOffsets[0] = 0; 4357 offsetsCopy[0] = 0; 4358 if (numFields) { 4359 PetscInt f; 4360 4361 for (f = 0; f < numFields; f++) { 4362 PetscInt fDof; 4363 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 4364 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4365 } 4366 ierr = DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 4367 } else { 4368 ierr = DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 4369 rowOffsets[1] = offsetsCopy[0]; 4370 } 4371 4372 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 4373 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 4374 leafEnd = leafStart + numLeaves; 4375 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4376 for (l = leafStart; l < leafEnd; l++) { 4377 PetscInt numIndices, childId, offset; 4378 const PetscScalar *childValues; 4379 4380 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 4381 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 4382 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4383 childValues = &rootValues[offset]; 4384 numIndices--; 4385 4386 if (childId == -2) { /* skip */ 4387 continue; 4388 } else if (childId == -1) { /* equivalent points: scatter */ 4389 PetscInt m; 4390 4391 contribute = PETSC_TRUE; 4392 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4393 } else { /* contributions from children: sum with injectors from reference tree */ 4394 PetscInt parentId, f, lim; 4395 4396 contribute = PETSC_TRUE; 4397 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 4398 4399 lim = PetscMax(1,numFields); 4400 offsets[0] = 0; 4401 if (numFields) { 4402 PetscInt f; 4403 4404 for (f = 0; f < numFields; f++) { 4405 PetscInt fDof; 4406 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4407 4408 offsets[f + 1] = fDof + offsets[f]; 4409 } 4410 } 4411 else { 4412 PetscInt cDof; 4413 4414 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4415 offsets[1] = cDof; 4416 } 4417 for (f = 0; f < lim; f++) { 4418 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4419 PetscInt n = offsets[f+1]-offsets[f]; 4420 PetscInt m = rowOffsets[f+1]-rowOffsets[f]; 4421 PetscInt i, j; 4422 const PetscScalar *colValues = &childValues[offsets[f]]; 4423 4424 for (i = 0; i < m; i++) { 4425 PetscScalar val = 0.; 4426 for (j = 0; j < n; j++) { 4427 val += childMat[n * i + j] * colValues[j]; 4428 } 4429 parentValues[rowOffsets[f] + i] += val; 4430 } 4431 } 4432 } 4433 } 4434 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);} 4435 } 4436 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4437 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4438 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr); 4439 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4440 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4441 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4442 PetscFunctionReturn(0); 4443 } 4444 4445 /*@ 4446 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4447 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4448 coarsening and refinement at the same time. 4449 4450 collective 4451 4452 Input Parameters: 4453 + dmIn - The DMPlex mesh for the input vector 4454 . vecIn - The input vector 4455 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4456 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4457 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4458 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4459 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4460 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4461 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4462 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4463 point j in dmOut is not a leaf of sfRefine. 4464 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4465 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4466 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4467 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4468 - time - Used if boundary values are time dependent. 4469 4470 Output Parameters: 4471 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered 4472 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4473 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4474 coarse points to fine points. 4475 4476 Level: developer 4477 4478 .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients() 4479 @*/ 4480 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4481 { 4482 PetscErrorCode ierr; 4483 4484 PetscFunctionBegin; 4485 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4486 if (sfRefine) { 4487 Vec vecInLocal; 4488 DM dmGrad = NULL; 4489 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4490 4491 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4492 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4493 { 4494 PetscInt numFields, i; 4495 4496 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4497 for (i = 0; i < numFields; i++) { 4498 PetscObject obj; 4499 PetscClassId classid; 4500 4501 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4502 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4503 if (classid == PETSCFV_CLASSID) { 4504 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4505 break; 4506 } 4507 } 4508 } 4509 if (useBCs) { 4510 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4511 } 4512 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4513 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4514 if (dmGrad) { 4515 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4516 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4517 } 4518 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr); 4519 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4520 if (dmGrad) { 4521 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4522 } 4523 } 4524 if (sfCoarsen) { 4525 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr); 4526 } 4527 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4528 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4529 PetscFunctionReturn(0); 4530 } 4531