1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petsc/private/petscfeimpl.h> 4 #include <petscsf.h> 5 #include <petscds.h> 6 7 /** hierarchy routines */ 8 9 /*@ 10 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 11 12 Not collective 13 14 Input Parameters: 15 + dm - The DMPlex object 16 - ref - The reference tree DMPlex object 17 18 Level: intermediate 19 20 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 21 @*/ 22 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 23 { 24 DM_Plex *mesh = (DM_Plex *)dm->data; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 if (ref) {PetscValidHeaderSpecific(ref, DM_CLASSID, 2);} 30 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 31 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 32 mesh->referenceTree = ref; 33 PetscFunctionReturn(0); 34 } 35 36 /*@ 37 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 38 39 Not collective 40 41 Input Parameters: 42 . dm - The DMPlex object 43 44 Output Parameters 45 . ref - The reference tree DMPlex object 46 47 Level: intermediate 48 49 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 50 @*/ 51 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 52 { 53 DM_Plex *mesh = (DM_Plex *)dm->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 57 PetscValidPointer(ref,2); 58 *ref = mesh->referenceTree; 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 63 { 64 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 if (parentOrientA == parentOrientB) { 69 if (childOrientB) *childOrientB = childOrientA; 70 if (childB) *childB = childA; 71 PetscFunctionReturn(0); 72 } 73 for (dim = 0; dim < 3; dim++) { 74 ierr = DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);CHKERRQ(ierr); 75 if (parent >= dStart && parent <= dEnd) { 76 break; 77 } 78 } 79 if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim); 80 if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children"); 81 if (childA < dStart || childA >= dEnd) { 82 /* this is a lower-dimensional child: bootstrap */ 83 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 84 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 85 86 ierr = DMPlexGetSupportSize(dm,childA,&size);CHKERRQ(ierr); 87 ierr = DMPlexGetSupport(dm,childA,&supp);CHKERRQ(ierr); 88 89 /* find a point sA in supp(childA) that has the same parent */ 90 for (i = 0; i < size; i++) { 91 PetscInt sParent; 92 93 sA = supp[i]; 94 if (sA == parent) continue; 95 ierr = DMPlexGetTreeParent(dm,sA,&sParent,NULL);CHKERRQ(ierr); 96 if (sParent == parent) { 97 break; 98 } 99 } 100 if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children"); 101 /* find out which point sB is in an equivalent position to sA under 102 * parentOrientB */ 103 ierr = DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);CHKERRQ(ierr); 104 ierr = DMPlexGetConeSize(dm,sA,&sConeSize);CHKERRQ(ierr); 105 ierr = DMPlexGetCone(dm,sA,&coneA);CHKERRQ(ierr); 106 ierr = DMPlexGetCone(dm,sB,&coneB);CHKERRQ(ierr); 107 ierr = DMPlexGetConeOrientation(dm,sA,&oA);CHKERRQ(ierr); 108 ierr = DMPlexGetConeOrientation(dm,sB,&oB);CHKERRQ(ierr); 109 /* step through the cone of sA in natural order */ 110 for (i = 0; i < sConeSize; i++) { 111 if (coneA[i] == childA) { 112 /* if childA is at position i in coneA, 113 * then we want the point that is at sOrientB*i in coneB */ 114 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize); 115 if (childB) *childB = coneB[j]; 116 if (childOrientB) { 117 PetscInt oBtrue; 118 119 ierr = DMPlexGetConeSize(dm,childA,&coneSize);CHKERRQ(ierr); 120 /* compose sOrientB and oB[j] */ 121 if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge"); 122 /* we may have to flip an edge */ 123 oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0; 124 ABswap = DihedralSwap(coneSize,oA[i],oBtrue);CHKERRQ(ierr); 125 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 126 } 127 break; 128 } 129 } 130 if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch"); 131 PetscFunctionReturn(0); 132 } 133 /* get the cone size and symmetry swap */ 134 ierr = DMPlexGetConeSize(dm,parent,&coneSize);CHKERRQ(ierr); 135 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 136 if (dim == 2) { 137 /* orientations refer to cones: we want them to refer to vertices: 138 * if it's a rotation, they are the same, but if the order is reversed, a 139 * permutation that puts side i first does *not* put vertex i first */ 140 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 141 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 142 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 143 } else { 144 ABswapVert = ABswap; 145 } 146 if (childB) { 147 /* assume that each child corresponds to a vertex, in the same order */ 148 PetscInt p, posA = -1, numChildren, i; 149 const PetscInt *children; 150 151 /* count which position the child is in */ 152 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 153 for (i = 0; i < numChildren; i++) { 154 p = children[i]; 155 if (p == childA) { 156 posA = i; 157 break; 158 } 159 } 160 if (posA >= coneSize) { 161 /* this is the triangle in the middle of a uniformly refined triangle: it 162 * is invariant */ 163 if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else"); 164 *childB = childA; 165 } 166 else { 167 /* figure out position B by applying ABswapVert */ 168 PetscInt posB; 169 170 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize); 171 if (childB) *childB = children[posB]; 172 } 173 } 174 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 175 PetscFunctionReturn(0); 176 } 177 178 /*@ 179 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another 180 181 Input Parameters: 182 + dm - the reference tree DMPlex object 183 . parent - the parent point 184 . parentOrientA - the reference orientation for describing the parent 185 . childOrientA - the reference orientation for describing the child 186 . childA - the reference childID for describing the child 187 - parentOrientB - the new orientation for describing the parent 188 189 Output Parameters: 190 + childOrientB - if not NULL, set to the new oreintation for describing the child 191 - childB - if not NULL, the new childID for describing the child 192 193 Level: developer 194 195 .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree() 196 @*/ 197 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 198 { 199 DM_Plex *mesh = (DM_Plex *)dm->data; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 204 if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented"); 205 ierr = mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool); 210 211 PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref) 212 { 213 MPI_Comm comm; 214 PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 215 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 216 DMLabel identity, identityRef; 217 PetscSection unionSection, unionConeSection, parentSection; 218 PetscScalar *unionCoords; 219 IS perm; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 comm = PetscObjectComm((PetscObject)K); 224 ierr = DMGetDimension(K, &dim);CHKERRQ(ierr); 225 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 226 ierr = DMGetLabel(K, labelName, &identity);CHKERRQ(ierr); 227 ierr = DMGetLabel(Kref, labelName, &identityRef);CHKERRQ(ierr); 228 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 229 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 230 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 231 /* count points that will go in the union */ 232 for (p = pStart; p < pEnd; p++) { 233 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 234 } 235 for (p = pRefStart; p < pRefEnd; p++) { 236 PetscInt q, qSize; 237 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 238 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 239 if (qSize > 1) { 240 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 241 } 242 } 243 ierr = PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);CHKERRQ(ierr); 244 offset = 0; 245 /* stratify points in the union by topological dimension */ 246 for (d = 0; d <= dim; d++) { 247 PetscInt cStart, cEnd, c; 248 249 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 250 for (c = cStart; c < cEnd; c++) { 251 permvals[offset++] = c; 252 } 253 254 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 255 for (c = cStart; c < cEnd; c++) { 256 permvals[offset++] = c + (pEnd - pStart); 257 } 258 } 259 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 260 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 261 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 262 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 263 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 264 /* count dimension points */ 265 for (d = 0; d <= dim; d++) { 266 PetscInt cStart, cOff, cOff2; 267 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 268 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 269 if (d < dim) { 270 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 271 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 272 } 273 else { 274 cOff2 = numUnionPoints; 275 } 276 numDimPoints[dim - d] = cOff2 - cOff; 277 } 278 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 279 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 280 /* count the cones in the union */ 281 for (p = pStart; p < pEnd; p++) { 282 PetscInt dof, uOff; 283 284 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 285 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 286 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 287 coneSizes[uOff] = dof; 288 } 289 for (p = pRefStart; p < pRefEnd; p++) { 290 PetscInt dof, uDof, uOff; 291 292 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 293 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 294 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 295 if (uDof) { 296 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 297 coneSizes[uOff] = dof; 298 } 299 } 300 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 301 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 302 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 303 /* write the cones in the union */ 304 for (p = pStart; p < pEnd; p++) { 305 PetscInt dof, uOff, c, cOff; 306 const PetscInt *cone, *orientation; 307 308 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 309 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 310 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 311 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 312 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 313 for (c = 0; c < dof; c++) { 314 PetscInt e, eOff; 315 e = cone[c]; 316 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 317 unionCones[cOff + c] = eOff; 318 unionOrientations[cOff + c] = orientation[c]; 319 } 320 } 321 for (p = pRefStart; p < pRefEnd; p++) { 322 PetscInt dof, uDof, uOff, c, cOff; 323 const PetscInt *cone, *orientation; 324 325 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 326 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 327 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 328 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 329 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 330 if (uDof) { 331 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 332 for (c = 0; c < dof; c++) { 333 PetscInt e, eOff, eDof; 334 335 e = cone[c]; 336 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 337 if (eDof) { 338 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 339 } 340 else { 341 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 342 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 343 } 344 unionCones[cOff + c] = eOff; 345 unionOrientations[cOff + c] = orientation[c]; 346 } 347 } 348 } 349 /* get the coordinates */ 350 { 351 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 352 PetscSection KcoordsSec, KrefCoordsSec; 353 Vec KcoordsVec, KrefCoordsVec; 354 PetscScalar *Kcoords; 355 356 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 357 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 358 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 359 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 360 361 numVerts = numDimPoints[0]; 362 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 363 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 364 365 offset = 0; 366 for (v = vStart; v < vEnd; v++) { 367 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 368 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 369 for (d = 0; d < dim; d++) { 370 unionCoords[offset * dim + d] = Kcoords[d]; 371 } 372 offset++; 373 } 374 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 375 for (v = vRefStart; v < vRefEnd; v++) { 376 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 377 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 378 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 379 if (vDof) { 380 for (d = 0; d < dim; d++) { 381 unionCoords[offset * dim + d] = Kcoords[d]; 382 } 383 offset++; 384 } 385 } 386 } 387 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 388 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 389 ierr = DMSetDimension(*ref,dim);CHKERRQ(ierr); 390 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 391 /* set the tree */ 392 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 393 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 394 for (p = pRefStart; p < pRefEnd; p++) { 395 PetscInt uDof, uOff; 396 397 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 398 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 399 if (uDof) { 400 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 401 } 402 } 403 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 404 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 405 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 406 for (p = pRefStart; p < pRefEnd; p++) { 407 PetscInt uDof, uOff; 408 409 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 410 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 411 if (uDof) { 412 PetscInt pOff, parent, parentU; 413 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 414 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 415 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 416 parents[pOff] = parentU; 417 childIDs[pOff] = uOff; 418 } 419 } 420 ierr = DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 421 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 422 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 423 424 /* clean up */ 425 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 426 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 427 ierr = ISDestroy(&perm);CHKERRQ(ierr); 428 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 429 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 430 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 431 PetscFunctionReturn(0); 432 } 433 434 /*@ 435 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 436 437 Collective on comm 438 439 Input Parameters: 440 + comm - the MPI communicator 441 . dim - the spatial dimension 442 - simplex - Flag for simplex, otherwise use a tensor-product cell 443 444 Output Parameters: 445 . ref - the reference tree DMPlex object 446 447 Level: intermediate 448 449 .keywords: reference cell 450 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 451 @*/ 452 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 453 { 454 DM_Plex *mesh; 455 DM K, Kref; 456 PetscInt p, pStart, pEnd; 457 DMLabel identity; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 #if 1 462 comm = PETSC_COMM_SELF; 463 #endif 464 /* create a reference element */ 465 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 466 ierr = DMCreateLabel(K, "identity");CHKERRQ(ierr); 467 ierr = DMGetLabel(K, "identity", &identity);CHKERRQ(ierr); 468 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 469 for (p = pStart; p < pEnd; p++) { 470 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 471 } 472 /* refine it */ 473 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 474 475 /* the reference tree is the union of these two, without duplicating 476 * points that appear in both */ 477 ierr = DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);CHKERRQ(ierr); 478 mesh = (DM_Plex *) (*ref)->data; 479 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default; 480 ierr = DMDestroy(&K);CHKERRQ(ierr); 481 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484 485 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 486 { 487 DM_Plex *mesh = (DM_Plex *)dm->data; 488 PetscSection childSec, pSec; 489 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 490 PetscInt *offsets, *children, pStart, pEnd; 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 495 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 496 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 497 pSec = mesh->parentSection; 498 if (!pSec) PetscFunctionReturn(0); 499 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 500 for (p = 0; p < pSize; p++) { 501 PetscInt par = mesh->parents[p]; 502 503 parMax = PetscMax(parMax,par+1); 504 parMin = PetscMin(parMin,par); 505 } 506 if (parMin > parMax) { 507 parMin = -1; 508 parMax = -1; 509 } 510 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 511 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 512 for (p = 0; p < pSize; p++) { 513 PetscInt par = mesh->parents[p]; 514 515 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 516 } 517 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 518 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 519 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 520 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 521 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 522 for (p = pStart; p < pEnd; p++) { 523 PetscInt dof, off, i; 524 525 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 526 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 527 for (i = 0; i < dof; i++) { 528 PetscInt par = mesh->parents[off + i], cOff; 529 530 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 531 children[cOff + offsets[par-parMin]++] = p; 532 } 533 } 534 mesh->childSection = childSec; 535 mesh->children = children; 536 ierr = PetscFree(offsets);CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 541 { 542 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 543 const PetscInt *vals; 544 PetscSection secNew; 545 PetscBool anyNew, globalAnyNew; 546 PetscBool compress; 547 PetscErrorCode ierr; 548 549 PetscFunctionBegin; 550 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 551 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 552 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 553 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 554 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 555 for (i = 0; i < size; i++) { 556 PetscInt dof; 557 558 p = vals[i]; 559 if (p < pStart || p >= pEnd) continue; 560 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 561 if (dof) break; 562 } 563 if (i == size) { 564 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 565 anyNew = PETSC_FALSE; 566 compress = PETSC_FALSE; 567 sizeNew = 0; 568 } 569 else { 570 anyNew = PETSC_TRUE; 571 for (p = pStart; p < pEnd; p++) { 572 PetscInt dof, off; 573 574 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 575 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 576 for (i = 0; i < dof; i++) { 577 PetscInt q = vals[off + i], qDof = 0; 578 579 if (q >= pStart && q < pEnd) { 580 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 581 } 582 if (qDof) { 583 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 584 } 585 else { 586 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 587 } 588 } 589 } 590 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 591 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 592 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 593 compress = PETSC_FALSE; 594 for (p = pStart; p < pEnd; p++) { 595 PetscInt dof, off, count, offNew, dofNew; 596 597 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 598 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 599 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 600 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 601 count = 0; 602 for (i = 0; i < dof; i++) { 603 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 604 605 if (q >= pStart && q < pEnd) { 606 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 607 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 608 } 609 if (qDof) { 610 PetscInt oldCount = count; 611 612 for (j = 0; j < qDof; j++) { 613 PetscInt k, r = vals[qOff + j]; 614 615 for (k = 0; k < oldCount; k++) { 616 if (valsNew[offNew + k] == r) { 617 break; 618 } 619 } 620 if (k == oldCount) { 621 valsNew[offNew + count++] = r; 622 } 623 } 624 } 625 else { 626 PetscInt k, oldCount = count; 627 628 for (k = 0; k < oldCount; k++) { 629 if (valsNew[offNew + k] == q) { 630 break; 631 } 632 } 633 if (k == oldCount) { 634 valsNew[offNew + count++] = q; 635 } 636 } 637 } 638 if (count < dofNew) { 639 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 640 compress = PETSC_TRUE; 641 } 642 } 643 } 644 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 645 ierr = MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 646 if (!globalAnyNew) { 647 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 648 *sectionNew = NULL; 649 *isNew = NULL; 650 } 651 else { 652 PetscBool globalCompress; 653 654 ierr = MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 655 if (compress) { 656 PetscSection secComp; 657 PetscInt *valsComp = NULL; 658 659 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 660 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 661 for (p = pStart; p < pEnd; p++) { 662 PetscInt dof; 663 664 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 665 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 666 } 667 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 668 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 669 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 670 for (p = pStart; p < pEnd; p++) { 671 PetscInt dof, off, offNew, j; 672 673 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 674 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 675 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 676 for (j = 0; j < dof; j++) { 677 valsComp[offNew + j] = valsNew[off + j]; 678 } 679 } 680 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 681 secNew = secComp; 682 ierr = PetscFree(valsNew);CHKERRQ(ierr); 683 valsNew = valsComp; 684 } 685 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 686 } 687 PetscFunctionReturn(0); 688 } 689 690 static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm) 691 { 692 PetscInt p, pStart, pEnd, *anchors, size; 693 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 694 PetscSection aSec; 695 DMLabel canonLabel; 696 IS aIS; 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 701 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 702 ierr = DMGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr); 703 for (p = pStart; p < pEnd; p++) { 704 PetscInt parent; 705 706 if (canonLabel) { 707 PetscInt canon; 708 709 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 710 if (p != canon) continue; 711 } 712 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 713 if (parent != p) { 714 aMin = PetscMin(aMin,p); 715 aMax = PetscMax(aMax,p+1); 716 } 717 } 718 if (aMin > aMax) { 719 aMin = -1; 720 aMax = -1; 721 } 722 ierr = PetscSectionCreate(PETSC_COMM_SELF,&aSec);CHKERRQ(ierr); 723 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 724 for (p = aMin; p < aMax; p++) { 725 PetscInt parent, ancestor = p; 726 727 if (canonLabel) { 728 PetscInt canon; 729 730 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 731 if (p != canon) continue; 732 } 733 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 734 while (parent != ancestor) { 735 ancestor = parent; 736 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 737 } 738 if (ancestor != p) { 739 PetscInt closureSize, *closure = NULL; 740 741 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 742 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 743 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 744 } 745 } 746 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 747 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 748 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 749 for (p = aMin; p < aMax; p++) { 750 PetscInt parent, ancestor = p; 751 752 if (canonLabel) { 753 PetscInt canon; 754 755 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 756 if (p != canon) continue; 757 } 758 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 759 while (parent != ancestor) { 760 ancestor = parent; 761 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 762 } 763 if (ancestor != p) { 764 PetscInt j, closureSize, *closure = NULL, aOff; 765 766 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 767 768 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 769 for (j = 0; j < closureSize; j++) { 770 anchors[aOff + j] = closure[2*j]; 771 } 772 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 773 } 774 } 775 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 776 { 777 PetscSection aSecNew = aSec; 778 IS aISNew = aIS; 779 780 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 781 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 782 while (aSecNew) { 783 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 784 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 785 aSec = aSecNew; 786 aIS = aISNew; 787 aSecNew = NULL; 788 aISNew = NULL; 789 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 790 } 791 } 792 ierr = DMPlexSetAnchors(dm,aSec,aIS);CHKERRQ(ierr); 793 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 794 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 798 static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp) 799 { 800 PetscErrorCode ierr; 801 802 PetscFunctionBegin; 803 if (numTrueSupp[p] == -1) { 804 PetscInt i, alldof; 805 const PetscInt *supp; 806 PetscInt count = 0; 807 808 ierr = DMPlexGetSupportSize(dm,p,&alldof);CHKERRQ(ierr); 809 ierr = DMPlexGetSupport(dm,p,&supp);CHKERRQ(ierr); 810 for (i = 0; i < alldof; i++) { 811 PetscInt q = supp[i], numCones, j; 812 const PetscInt *cone; 813 814 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 815 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 816 for (j = 0; j < numCones; j++) { 817 if (cone[j] == p) break; 818 } 819 if (j < numCones) count++; 820 } 821 numTrueSupp[p] = count; 822 } 823 *dof = numTrueSupp[p]; 824 PetscFunctionReturn(0); 825 } 826 827 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm) 828 { 829 DM_Plex *mesh = (DM_Plex *)dm->data; 830 PetscSection newSupportSection; 831 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth; 832 PetscInt *numTrueSupp; 833 PetscInt *offsets; 834 PetscErrorCode ierr; 835 836 PetscFunctionBegin; 837 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 838 /* symmetrize the hierarchy */ 839 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 840 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);CHKERRQ(ierr); 841 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 842 ierr = PetscSectionSetChart(newSupportSection,pStart,pEnd);CHKERRQ(ierr); 843 ierr = PetscCalloc1(pEnd,&offsets);CHKERRQ(ierr); 844 ierr = PetscMalloc1(pEnd,&numTrueSupp);CHKERRQ(ierr); 845 for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1; 846 /* if a point is in the (true) support of q, it should be in the support of 847 * parent(q) */ 848 for (d = 0; d <= depth; d++) { 849 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 850 for (p = pStart; p < pEnd; ++p) { 851 PetscInt dof, q, qdof, parent; 852 853 ierr = DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);CHKERRQ(ierr); 854 ierr = PetscSectionAddDof(newSupportSection, p, dof);CHKERRQ(ierr); 855 q = p; 856 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 857 while (parent != q && parent >= pStart && parent < pEnd) { 858 q = parent; 859 860 ierr = DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);CHKERRQ(ierr); 861 ierr = PetscSectionAddDof(newSupportSection,p,qdof);CHKERRQ(ierr); 862 ierr = PetscSectionAddDof(newSupportSection,q,dof);CHKERRQ(ierr); 863 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 864 } 865 } 866 } 867 ierr = PetscSectionSetUp(newSupportSection);CHKERRQ(ierr); 868 ierr = PetscSectionGetStorageSize(newSupportSection,&newSize);CHKERRQ(ierr); 869 ierr = PetscMalloc1(newSize,&newSupports);CHKERRQ(ierr); 870 for (d = 0; d <= depth; d++) { 871 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 872 for (p = pStart; p < pEnd; p++) { 873 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent; 874 875 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 876 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 877 ierr = PetscSectionGetDof(newSupportSection, p, &newDof);CHKERRQ(ierr); 878 ierr = PetscSectionGetOffset(newSupportSection, p, &newOff);CHKERRQ(ierr); 879 for (i = 0; i < dof; i++) { 880 PetscInt numCones, j; 881 const PetscInt *cone; 882 PetscInt q = mesh->supports[off + i]; 883 884 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 885 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 886 for (j = 0; j < numCones; j++) { 887 if (cone[j] == p) break; 888 } 889 if (j < numCones) newSupports[newOff+offsets[p]++] = q; 890 } 891 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof); 892 893 q = p; 894 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 895 while (parent != q && parent >= pStart && parent < pEnd) { 896 q = parent; 897 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr); 898 ierr = PetscSectionGetOffset(mesh->supportSection, q, &qoff);CHKERRQ(ierr); 899 ierr = PetscSectionGetOffset(newSupportSection, q, &newqOff);CHKERRQ(ierr); 900 for (i = 0; i < qdof; i++) { 901 PetscInt numCones, j; 902 const PetscInt *cone; 903 PetscInt r = mesh->supports[qoff + i]; 904 905 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 906 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 907 for (j = 0; j < numCones; j++) { 908 if (cone[j] == q) break; 909 } 910 if (j < numCones) newSupports[newOff+offsets[p]++] = r; 911 } 912 for (i = 0; i < dof; i++) { 913 PetscInt numCones, j; 914 const PetscInt *cone; 915 PetscInt r = mesh->supports[off + i]; 916 917 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 918 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 919 for (j = 0; j < numCones; j++) { 920 if (cone[j] == p) break; 921 } 922 if (j < numCones) newSupports[newqOff+offsets[q]++] = r; 923 } 924 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 925 } 926 } 927 } 928 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr); 929 mesh->supportSection = newSupportSection; 930 ierr = PetscFree(mesh->supports);CHKERRQ(ierr); 931 mesh->supports = newSupports; 932 ierr = PetscFree(offsets);CHKERRQ(ierr); 933 ierr = PetscFree(numTrueSupp);CHKERRQ(ierr); 934 935 PetscFunctionReturn(0); 936 } 937 938 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat); 939 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat); 940 941 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports) 942 { 943 DM_Plex *mesh = (DM_Plex *)dm->data; 944 DM refTree; 945 PetscInt size; 946 PetscErrorCode ierr; 947 948 PetscFunctionBegin; 949 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 950 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 951 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 952 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 953 mesh->parentSection = parentSection; 954 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 955 if (parents != mesh->parents) { 956 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 957 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 958 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 959 } 960 if (childIDs != mesh->childIDs) { 961 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 962 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 963 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 964 } 965 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 966 if (refTree) { 967 DMLabel canonLabel; 968 969 ierr = DMGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr); 970 if (canonLabel) { 971 PetscInt i; 972 973 for (i = 0; i < size; i++) { 974 PetscInt canon; 975 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr); 976 if (canon >= 0) { 977 mesh->childIDs[i] = canon; 978 } 979 } 980 } 981 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference; 982 } 983 else { 984 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct; 985 } 986 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 987 if (computeCanonical) { 988 PetscInt d, dim; 989 990 /* add the canonical label */ 991 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 992 ierr = DMCreateLabel(dm,"canonical");CHKERRQ(ierr); 993 for (d = 0; d <= dim; d++) { 994 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 995 const PetscInt *cChildren; 996 997 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 998 for (p = dStart; p < dEnd; p++) { 999 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr); 1000 if (cNumChildren) { 1001 canon = p; 1002 break; 1003 } 1004 } 1005 if (canon == -1) continue; 1006 for (p = dStart; p < dEnd; p++) { 1007 PetscInt numChildren, i; 1008 const PetscInt *children; 1009 1010 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 1011 if (numChildren) { 1012 if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren); 1013 ierr = DMSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr); 1014 for (i = 0; i < numChildren; i++) { 1015 ierr = DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr); 1016 } 1017 } 1018 } 1019 } 1020 } 1021 if (exchangeSupports) { 1022 ierr = DMPlexTreeExchangeSupports(dm);CHKERRQ(ierr); 1023 } 1024 mesh->createanchors = DMPlexCreateAnchors_Tree; 1025 /* reset anchors */ 1026 ierr = DMPlexSetAnchors(dm,NULL,NULL);CHKERRQ(ierr); 1027 PetscFunctionReturn(0); 1028 } 1029 1030 /*@ 1031 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 1032 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 1033 tree root. 1034 1035 Collective on dm 1036 1037 Input Parameters: 1038 + dm - the DMPlex object 1039 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1040 offset indexes the parent and childID list; the reference count of parentSection is incremented 1041 . parents - a list of the point parents; copied, can be destroyed 1042 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1043 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 1044 1045 Level: intermediate 1046 1047 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1048 @*/ 1049 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 1050 { 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 /*@ 1059 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 1060 Collective on dm 1061 1062 Input Parameters: 1063 . dm - the DMPlex object 1064 1065 Output Parameters: 1066 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1067 offset indexes the parent and childID list 1068 . parents - a list of the point parents 1069 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1070 the child corresponds to the point in the reference tree with index childID 1071 . childSection - the inverse of the parent section 1072 - children - a list of the point children 1073 1074 Level: intermediate 1075 1076 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1077 @*/ 1078 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 1079 { 1080 DM_Plex *mesh = (DM_Plex *)dm->data; 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1084 if (parentSection) *parentSection = mesh->parentSection; 1085 if (parents) *parents = mesh->parents; 1086 if (childIDs) *childIDs = mesh->childIDs; 1087 if (childSection) *childSection = mesh->childSection; 1088 if (children) *children = mesh->children; 1089 PetscFunctionReturn(0); 1090 } 1091 1092 /*@ 1093 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve 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 Sieve 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 nFunctionals, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints) 1187 { 1188 PetscInt i, j, k, offset, qPoints; 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 ierr = PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);CHKERRQ(ierr); 1193 for (i = 0, offset = 0; i < nFunctionals; i++) { 1194 qPoints = pointsPerFn[i]; 1195 for (j = 0; j < nFunctionals; j++) { 1196 PetscScalar val = 0.; 1197 1198 for (k = 0; k < qPoints; k++) { 1199 val += work[(offset + k) * nFunctionals + j] * weights[k]; 1200 } 1201 ierr = MatSetValue(basisAtPoints,j,i,val,INSERT_VALUES);CHKERRQ(ierr); 1202 } 1203 offset += qPoints; 1204 } 1205 ierr = MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1206 ierr = MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat) 1211 { 1212 PetscDS ds; 1213 PetscInt spdim; 1214 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 1215 const PetscInt *anchors; 1216 PetscSection aSec; 1217 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 1218 IS aIS; 1219 PetscErrorCode ierr; 1220 1221 PetscFunctionBegin; 1222 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1223 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1224 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1225 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1226 ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr); 1227 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 1228 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 1229 ierr = DMGetDimension(dm,&spdim);CHKERRQ(ierr); 1230 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 1231 1232 for (f = 0; f < numFields; f++) { 1233 PetscObject disc; 1234 PetscClassId id; 1235 PetscSpace bspace; 1236 PetscDualSpace dspace; 1237 PetscInt i, j, k, nPoints, offset; 1238 PetscInt fSize, fComp, maxDof; 1239 PetscReal *weights, *pointsRef, *pointsReal, *work; 1240 PetscScalar *scwork, *X; 1241 PetscInt *sizes, *workIndRow, *workIndCol; 1242 Mat Amat, Bmat, Xmat; 1243 const PetscInt *numDof = NULL; 1244 const PetscInt ***perms = NULL; 1245 const PetscScalar ***flips = NULL; 1246 1247 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 1248 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1249 if (id == PETSCFE_CLASSID) { 1250 PetscFE fe = (PetscFE) disc; 1251 1252 ierr = PetscFEGetBasisSpace(fe,&bspace);CHKERRQ(ierr); 1253 ierr = PetscFEGetDualSpace(fe,&dspace);CHKERRQ(ierr); 1254 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr); 1255 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1256 } 1257 else if (id == PETSCFV_CLASSID) { 1258 PetscFV fv = (PetscFV) disc; 1259 1260 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);CHKERRQ(ierr); 1261 ierr = PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1262 ierr = PetscSpaceSetOrder(bspace,0);CHKERRQ(ierr); 1263 ierr = PetscSpacePolynomialSetNumVariables(bspace,spdim);CHKERRQ(ierr); 1264 ierr = PetscSpaceSetUp(bspace);CHKERRQ(ierr); 1265 ierr = PetscFVGetDualSpace(fv,&dspace);CHKERRQ(ierr); 1266 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr); 1267 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 1268 } 1269 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 1270 ierr = PetscDualSpaceGetNumDof(dspace,&numDof);CHKERRQ(ierr); 1271 for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);} 1272 ierr = PetscDualSpaceGetSymmetries(dspace,&perms,&flips);CHKERRQ(ierr); 1273 1274 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 1275 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 1276 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 1277 ierr = MatSetUp(Amat);CHKERRQ(ierr); 1278 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 1279 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 1280 nPoints = 0; 1281 for (i = 0; i < fSize; i++) { 1282 PetscInt qPoints; 1283 PetscQuadrature quad; 1284 1285 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr); 1286 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1287 nPoints += qPoints; 1288 } 1289 ierr = PetscMalloc7(fSize,&sizes,nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,fSize*nPoints,&work,maxDof,&workIndRow,maxDof,&workIndCol);CHKERRQ(ierr); 1290 ierr = PetscMalloc1(maxDof * maxDof,&scwork);CHKERRQ(ierr); 1291 offset = 0; 1292 for (i = 0; i < fSize; i++) { 1293 PetscInt qPoints; 1294 const PetscReal *p, *w; 1295 PetscQuadrature quad; 1296 1297 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr); 1298 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 1299 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 1300 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 1301 sizes[i] = qPoints; 1302 offset += qPoints; 1303 } 1304 ierr = EvaluateBasis(bspace,fSize,nPoints,sizes,pointsRef,weights,work,Amat);CHKERRQ(ierr); 1305 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1306 for (c = cStart; c < cEnd; c++) { 1307 PetscInt parent; 1308 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 1309 PetscInt *childOffsets, *parentOffsets; 1310 1311 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 1312 if (parent == c) continue; 1313 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1314 for (i = 0; i < closureSize; i++) { 1315 PetscInt p = closure[2*i]; 1316 PetscInt conDof; 1317 1318 if (p < conStart || p >= conEnd) continue; 1319 if (numFields) { 1320 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1321 } 1322 else { 1323 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1324 } 1325 if (conDof) break; 1326 } 1327 if (i == closureSize) { 1328 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1329 continue; 1330 } 1331 1332 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 1333 ierr = DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 1334 for (i = 0; i < nPoints; i++) { 1335 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp); 1336 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]); 1337 } 1338 ierr = EvaluateBasis(bspace,fSize,nPoints,sizes,pointsReal,weights,work,Bmat);CHKERRQ(ierr); 1339 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1340 ierr = MatDenseGetArray(Xmat,&X);CHKERRQ(ierr); 1341 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1342 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1343 childOffsets[0] = 0; 1344 for (i = 0; i < closureSize; i++) { 1345 PetscInt p = closure[2*i]; 1346 PetscInt dof; 1347 1348 if (numFields) { 1349 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1350 } 1351 else { 1352 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1353 } 1354 childOffsets[i+1]=childOffsets[i]+dof / fComp; 1355 } 1356 parentOffsets[0] = 0; 1357 for (i = 0; i < closureSizeP; i++) { 1358 PetscInt p = closureP[2*i]; 1359 PetscInt dof; 1360 1361 if (numFields) { 1362 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1363 } 1364 else { 1365 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1366 } 1367 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 1368 } 1369 for (i = 0; i < closureSize; i++) { 1370 PetscInt conDof, conOff, aDof, aOff, nWork; 1371 PetscInt p = closure[2*i]; 1372 PetscInt o = closure[2*i+1]; 1373 const PetscInt *perm; 1374 const PetscScalar *flip; 1375 1376 if (p < conStart || p >= conEnd) continue; 1377 if (numFields) { 1378 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1379 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1380 } 1381 else { 1382 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1383 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1384 } 1385 if (!conDof) continue; 1386 perm = (perms && perms[i]) ? perms[i][o] : NULL; 1387 flip = (flips && flips[i]) ? flips[i][o] : NULL; 1388 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1389 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1390 nWork = childOffsets[i+1]-childOffsets[i]; 1391 for (k = 0; k < aDof; k++) { 1392 PetscInt a = anchors[aOff + k]; 1393 PetscInt aSecDof, aSecOff; 1394 1395 if (numFields) { 1396 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1397 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1398 } 1399 else { 1400 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1401 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1402 } 1403 if (!aSecDof) continue; 1404 1405 for (j = 0; j < closureSizeP; j++) { 1406 PetscInt q = closureP[2*j]; 1407 PetscInt oq = closureP[2*j+1]; 1408 1409 if (q == a) { 1410 PetscInt r, s, t, nWorkP; 1411 const PetscInt *permP; 1412 const PetscScalar *flipP; 1413 1414 permP = (perms && perms[j]) ? perms[j][oq] : NULL; 1415 flipP = (flips && flips[j]) ? flips[j][oq] : NULL; 1416 nWorkP = parentOffsets[j+1]-parentOffsets[j]; 1417 /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the 1418 * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArray is 1419 * column-major, so transpose-transpose = do nothing */ 1420 for (r = 0; r < nWork; r++) { 1421 for (s = 0; s < nWorkP; s++) { 1422 scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])]; 1423 } 1424 } 1425 for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + fComp * r;} 1426 for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + fComp * s;} 1427 if (flip) { 1428 for (r = 0; r < nWork; r++) { 1429 for (s = 0; s < nWorkP; s++) { 1430 scwork[r * nWorkP + s] *= flip[r]; 1431 } 1432 } 1433 } 1434 if (flipP) { 1435 for (r = 0; r < nWork; r++) { 1436 for (s = 0; s < nWorkP; s++) { 1437 scwork[r * nWorkP + s] *= flipP[s]; 1438 } 1439 } 1440 } 1441 for (t = 0; t < fComp; t++) { 1442 ierr = MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);CHKERRQ(ierr); 1443 if (t < fComp - 1) { 1444 for (r = 0; r < nWork; r++) {workIndRow[r]++;} 1445 for (s = 0; s < nWorkP; s++) {workIndCol[s]++;} 1446 } 1447 } 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 = DMGetDefaultSection(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);CHKERRQ(ierr); 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 2112 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr); 2113 if (kParent != v) { 2114 /* this is a new vertex */ 2115 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr); 2116 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2117 CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);CHKERRQ(ierr); 2118 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2119 offset += dim; 2120 } 2121 } 2122 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2123 } 2124 2125 /* need to reverse the order of pNewCount: vertices first, cells last */ 2126 for (d = 0; d < (dim + 1) / 2; d++) { 2127 PetscInt tmp; 2128 2129 tmp = pNewCount[d]; 2130 pNewCount[d] = pNewCount[dim - d]; 2131 pNewCount[dim - d] = tmp; 2132 } 2133 2134 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 2135 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2136 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 2137 2138 /* clean up */ 2139 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 2140 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 2141 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 2142 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 2143 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 2144 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 2145 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 2146 } 2147 else { 2148 PetscInt p, counts[4]; 2149 PetscInt *coneSizes, *cones, *orientations; 2150 Vec coordVec; 2151 PetscScalar *coords; 2152 2153 for (d = 0; d <= dim; d++) { 2154 PetscInt dStart, dEnd; 2155 2156 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 2157 counts[d] = dEnd - dStart; 2158 } 2159 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 2160 for (p = pStart; p < pEnd; p++) { 2161 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 2162 } 2163 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2164 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2165 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2166 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2167 2168 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2169 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2170 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2171 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2172 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2173 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2174 } 2175 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2176 2177 PetscFunctionReturn(0); 2178 } 2179 2180 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2181 { 2182 PetscSF coarseToFineEmbedded; 2183 PetscSection globalCoarse, globalFine; 2184 PetscSection localCoarse, localFine; 2185 PetscSection aSec, cSec; 2186 PetscSection rootIndicesSec, rootMatricesSec; 2187 PetscSection leafIndicesSec, leafMatricesSec; 2188 PetscInt *rootIndices, *leafIndices; 2189 PetscScalar *rootMatrices, *leafMatrices; 2190 IS aIS; 2191 const PetscInt *anchors; 2192 Mat cMat; 2193 PetscInt numFields, maxFields; 2194 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2195 PetscInt aStart, aEnd, cStart, cEnd; 2196 PetscInt *maxChildIds; 2197 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2198 const PetscInt ***perms; 2199 const PetscScalar ***flips; 2200 PetscErrorCode ierr; 2201 2202 PetscFunctionBegin; 2203 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 2204 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 2205 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 2206 { /* winnow fine points that don't have global dofs out of the sf */ 2207 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2208 const PetscInt *leaves; 2209 2210 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 2211 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2212 p = leaves ? leaves[l] : l; 2213 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2214 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2215 if ((dof - cdof) > 0) { 2216 numPointsWithDofs++; 2217 } 2218 } 2219 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 2220 for (l = 0, offset = 0; l < nleaves; l++) { 2221 p = leaves ? leaves[l] : l; 2222 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2223 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2224 if ((dof - cdof) > 0) { 2225 pointsWithDofs[offset++] = l; 2226 } 2227 } 2228 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 2229 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 2230 } 2231 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2232 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 2233 for (p = pStartC; p < pEndC; p++) { 2234 maxChildIds[p - pStartC] = -2; 2235 } 2236 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2237 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2238 2239 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 2240 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 2241 2242 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 2243 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 2244 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 2245 2246 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 2247 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 2248 2249 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2250 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 2251 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr); 2252 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr); 2253 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr); 2254 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 2255 maxFields = PetscMax(1,numFields); 2256 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); 2257 ierr = PetscMalloc2(maxFields+1,&perms,maxFields+1,&flips);CHKERRQ(ierr); 2258 ierr = PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));CHKERRQ(ierr); 2259 ierr = PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));CHKERRQ(ierr); 2260 2261 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2262 PetscInt dof, matSize = 0; 2263 PetscInt aDof = 0; 2264 PetscInt cDof = 0; 2265 PetscInt maxChildId = maxChildIds[p - pStartC]; 2266 PetscInt numRowIndices = 0; 2267 PetscInt numColIndices = 0; 2268 PetscInt f; 2269 2270 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2271 if (dof < 0) { 2272 dof = -(dof + 1); 2273 } 2274 if (p >= aStart && p < aEnd) { 2275 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2276 } 2277 if (p >= cStart && p < cEnd) { 2278 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr); 2279 } 2280 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2281 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2282 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2283 PetscInt *closure = NULL, closureSize, cl; 2284 2285 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2286 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2287 PetscInt c = closure[2 * cl], clDof; 2288 2289 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 2290 numRowIndices += clDof; 2291 for (f = 0; f < numFields; f++) { 2292 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr); 2293 offsets[f + 1] += clDof; 2294 } 2295 } 2296 for (f = 0; f < numFields; f++) { 2297 offsets[f + 1] += offsets[f]; 2298 newOffsets[f + 1] = offsets[f + 1]; 2299 } 2300 /* get the number of indices needed and their field offsets */ 2301 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2302 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2303 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2304 numColIndices = numRowIndices; 2305 matSize = 0; 2306 } 2307 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2308 matSize = 0; 2309 for (f = 0; f < numFields; f++) { 2310 PetscInt numRow, numCol; 2311 2312 numRow = offsets[f + 1] - offsets[f]; 2313 numCol = newOffsets[f + 1] - newOffsets[f]; 2314 matSize += numRow * numCol; 2315 } 2316 } 2317 else { 2318 matSize = numRowIndices * numColIndices; 2319 } 2320 } else if (maxChildId == -1) { 2321 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2322 PetscInt aOff, a; 2323 2324 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2325 for (f = 0; f < numFields; f++) { 2326 PetscInt fDof; 2327 2328 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2329 offsets[f+1] = fDof; 2330 } 2331 for (a = 0; a < aDof; a++) { 2332 PetscInt anchor = anchors[a + aOff], aLocalDof; 2333 2334 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr); 2335 numColIndices += aLocalDof; 2336 for (f = 0; f < numFields; f++) { 2337 PetscInt fDof; 2338 2339 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2340 newOffsets[f+1] += fDof; 2341 } 2342 } 2343 if (numFields) { 2344 matSize = 0; 2345 for (f = 0; f < numFields; f++) { 2346 matSize += offsets[f+1] * newOffsets[f+1]; 2347 } 2348 } 2349 else { 2350 matSize = numColIndices * dof; 2351 } 2352 } 2353 else { /* no children, and no constraints on dofs: just get the global indices */ 2354 numColIndices = dof; 2355 matSize = 0; 2356 } 2357 } 2358 /* we will pack the column indices with the field offsets */ 2359 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr); 2360 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr); 2361 } 2362 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr); 2363 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr); 2364 { 2365 PetscInt numRootIndices, numRootMatrices; 2366 2367 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 2368 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr); 2369 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr); 2370 for (p = pStartC; p < pEndC; p++) { 2371 PetscInt numRowIndices, numColIndices, matSize, dof; 2372 PetscInt pIndOff, pMatOff, f; 2373 PetscInt *pInd; 2374 PetscInt maxChildId = maxChildIds[p - pStartC]; 2375 PetscScalar *pMat = NULL; 2376 2377 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2378 if (!numColIndices) { 2379 continue; 2380 } 2381 for (f = 0; f <= numFields; f++) { 2382 offsets[f] = 0; 2383 newOffsets[f] = 0; 2384 offsetsCopy[f] = 0; 2385 newOffsetsCopy[f] = 0; 2386 } 2387 numColIndices -= 2 * numFields; 2388 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2389 pInd = &(rootIndices[pIndOff]); 2390 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr); 2391 if (matSize) { 2392 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2393 pMat = &rootMatrices[pMatOff]; 2394 } 2395 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2396 if (dof < 0) { 2397 dof = -(dof + 1); 2398 } 2399 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2400 PetscInt i, j; 2401 PetscInt numRowIndices = matSize / numColIndices; 2402 2403 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2404 PetscInt numIndices, *indices; 2405 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2406 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2407 for (i = 0; i < numColIndices; i++) { 2408 pInd[i] = indices[i]; 2409 } 2410 for (i = 0; i < numFields; i++) { 2411 pInd[numColIndices + i] = offsets[i+1]; 2412 pInd[numColIndices + numFields + i] = offsets[i+1]; 2413 } 2414 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2415 } 2416 else { 2417 PetscInt closureSize, *closure = NULL, cl; 2418 PetscScalar *pMatIn, *pMatModified; 2419 PetscInt numPoints,*points; 2420 2421 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2422 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2423 for (j = 0; j < numRowIndices; j++) { 2424 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2425 } 2426 } 2427 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2428 for (f = 0; f < maxFields; f++) { 2429 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2430 else {ierr = PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2431 } 2432 if (numFields) { 2433 for (cl = 0; cl < closureSize; cl++) { 2434 PetscInt c = closure[2 * cl]; 2435 2436 for (f = 0; f < numFields; f++) { 2437 PetscInt fDof; 2438 2439 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr); 2440 offsets[f + 1] += fDof; 2441 } 2442 } 2443 for (f = 0; f < numFields; f++) { 2444 offsets[f + 1] += offsets[f]; 2445 newOffsets[f + 1] = offsets[f + 1]; 2446 } 2447 } 2448 /* TODO : flips here ? */ 2449 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2450 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2451 for (f = 0; f < maxFields; f++) { 2452 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2453 else {ierr = PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2454 } 2455 for (f = 0; f < maxFields; f++) { 2456 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2457 else {ierr = PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2458 } 2459 if (!numFields) { 2460 for (i = 0; i < numRowIndices * numColIndices; i++) { 2461 pMat[i] = pMatModified[i]; 2462 } 2463 } 2464 else { 2465 PetscInt i, j, count; 2466 for (f = 0, count = 0; f < numFields; f++) { 2467 for (i = offsets[f]; i < offsets[f+1]; i++) { 2468 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2469 pMat[count] = pMatModified[i * numColIndices + j]; 2470 } 2471 } 2472 } 2473 } 2474 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);CHKERRQ(ierr); 2475 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2476 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2477 if (numFields) { 2478 for (f = 0; f < numFields; f++) { 2479 pInd[numColIndices + f] = offsets[f+1]; 2480 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2481 } 2482 for (cl = 0; cl < numPoints; cl++) { 2483 PetscInt globalOff, c = points[2*cl]; 2484 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2485 DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd); 2486 } 2487 } else { 2488 for (cl = 0; cl < numPoints; cl++) { 2489 PetscInt c = points[2*cl], globalOff; 2490 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2491 2492 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2493 DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd); 2494 } 2495 } 2496 for (f = 0; f < maxFields; f++) { 2497 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2498 else {ierr = PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2499 } 2500 ierr = DMRestoreWorkArray(coarse,numPoints,PETSC_SCALAR,&points);CHKERRQ(ierr); 2501 } 2502 } 2503 else if (matSize) { 2504 PetscInt cOff; 2505 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2506 2507 numRowIndices = matSize / numColIndices; 2508 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2509 ierr = DMGetWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2510 ierr = DMGetWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2511 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr); 2512 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2513 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2514 if (numFields) { 2515 for (f = 0; f < numFields; f++) { 2516 PetscInt fDof; 2517 2518 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr); 2519 offsets[f + 1] = fDof; 2520 for (a = 0; a < aDof; a++) { 2521 PetscInt anchor = anchors[a + aOff]; 2522 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2523 newOffsets[f + 1] += fDof; 2524 } 2525 } 2526 for (f = 0; f < numFields; f++) { 2527 offsets[f + 1] += offsets[f]; 2528 offsetsCopy[f + 1] = offsets[f + 1]; 2529 newOffsets[f + 1] += newOffsets[f]; 2530 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2531 } 2532 DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);CHKERRQ(ierr); 2533 for (a = 0; a < aDof; a++) { 2534 PetscInt anchor = anchors[a + aOff], lOff; 2535 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2536 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);CHKERRQ(ierr); 2537 } 2538 } 2539 else { 2540 DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);CHKERRQ(ierr); 2541 for (a = 0; a < aDof; a++) { 2542 PetscInt anchor = anchors[a + aOff], lOff; 2543 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2544 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);CHKERRQ(ierr); 2545 } 2546 } 2547 if (numFields) { 2548 PetscInt count, a; 2549 2550 for (f = 0, count = 0; f < numFields; f++) { 2551 PetscInt iSize = offsets[f + 1] - offsets[f]; 2552 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2553 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2554 count += iSize * jSize; 2555 pInd[numColIndices + f] = offsets[f+1]; 2556 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2557 } 2558 for (a = 0; a < aDof; a++) { 2559 PetscInt anchor = anchors[a + aOff]; 2560 PetscInt gOff; 2561 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2562 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2563 } 2564 } 2565 else { 2566 PetscInt a; 2567 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2568 for (a = 0; a < aDof; a++) { 2569 PetscInt anchor = anchors[a + aOff]; 2570 PetscInt gOff; 2571 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2572 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2573 } 2574 } 2575 ierr = DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2576 ierr = DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2577 } 2578 else { 2579 PetscInt gOff; 2580 2581 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 2582 if (numFields) { 2583 for (f = 0; f < numFields; f++) { 2584 PetscInt fDof; 2585 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2586 offsets[f + 1] = fDof + offsets[f]; 2587 } 2588 for (f = 0; f < numFields; f++) { 2589 pInd[numColIndices + f] = offsets[f+1]; 2590 pInd[numColIndices + numFields + f] = offsets[f+1]; 2591 } 2592 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2593 } 2594 else { 2595 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 = DMGetDefaultSection(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,PETSC_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 2840 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2841 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2842 if ((gDof - gcDof) <= 0) { 2843 continue; 2844 } 2845 childId = childIds[p-pStartF]; 2846 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2847 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2848 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2849 numColIndices -= 2 * numFields; 2850 pInd = &leafIndices[pIndOff]; 2851 offsets[0] = 0; 2852 offsetsCopy[0] = 0; 2853 newOffsets[0] = 0; 2854 newOffsetsCopy[0] = 0; 2855 rowOffsets[0] = 0; 2856 if (numFields) { 2857 PetscInt f; 2858 for (f = 0; f < numFields; f++) { 2859 PetscInt rowDof; 2860 2861 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2862 offsets[f + 1] = offsets[f] + rowDof; 2863 offsetsCopy[f + 1] = offsets[f + 1]; 2864 rowOffsets[f + 1] = pInd[numColIndices + f]; 2865 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2866 } 2867 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2868 } 2869 else { 2870 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2871 } 2872 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2873 if (!matSize) { /* incoming matrix is identity */ 2874 if (childId < 0) { /* no child interpolation: scatter */ 2875 if (numFields) { 2876 PetscInt f; 2877 for (f = 0; f < numFields; f++) { 2878 PetscInt numRows = offsets[f+1] - offsets[f], row; 2879 for (row = 0; row < numRows; row++) { 2880 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2881 } 2882 } 2883 } 2884 else { 2885 PetscInt numRows = gDof, row; 2886 for (row = 0; row < numRows; row++) { 2887 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2888 } 2889 } 2890 } 2891 else { /* interpolate from all */ 2892 if (numFields) { 2893 PetscInt f; 2894 for (f = 0; f < numFields; f++) { 2895 PetscInt numRows = offsets[f+1] - offsets[f]; 2896 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2897 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2898 } 2899 } 2900 else { 2901 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2902 } 2903 } 2904 } 2905 else { /* interpolate from all */ 2906 PetscInt pMatOff; 2907 PetscScalar *pMat; 2908 2909 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2910 pMat = &leafMatrices[pMatOff]; 2911 if (childId < 0) { /* copy the incoming matrix */ 2912 if (numFields) { 2913 PetscInt f, count; 2914 for (f = 0, count = 0; f < numFields; f++) { 2915 PetscInt numRows = offsets[f+1]-offsets[f]; 2916 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2917 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2918 PetscScalar *inMat = &pMat[count]; 2919 2920 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2921 count += numCols * numInRows; 2922 } 2923 } 2924 else { 2925 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2926 } 2927 } 2928 else { /* multiply the incoming matrix by the child interpolation */ 2929 if (numFields) { 2930 PetscInt f, count; 2931 for (f = 0, count = 0; f < numFields; f++) { 2932 PetscInt numRows = offsets[f+1]-offsets[f]; 2933 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2934 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2935 PetscScalar *inMat = &pMat[count]; 2936 PetscInt i, j, k; 2937 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2938 for (i = 0; i < numRows; i++) { 2939 for (j = 0; j < numCols; j++) { 2940 PetscScalar val = 0.; 2941 for (k = 0; k < numInRows; k++) { 2942 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2943 } 2944 pointWork[i * numCols + j] = val; 2945 } 2946 } 2947 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2948 count += numCols * numInRows; 2949 } 2950 } 2951 else { /* every dof gets a full row */ 2952 PetscInt numRows = gDof; 2953 PetscInt numCols = numColIndices; 2954 PetscInt numInRows = matSize / numColIndices; 2955 PetscInt i, j, k; 2956 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2957 for (i = 0; i < numRows; i++) { 2958 for (j = 0; j < numCols; j++) { 2959 PetscScalar val = 0.; 2960 for (k = 0; k < numInRows; k++) { 2961 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2962 } 2963 pointWork[i * numCols + j] = val; 2964 } 2965 } 2966 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 2967 } 2968 } 2969 } 2970 } 2971 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2972 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2973 ierr = PetscFree(pointWork);CHKERRQ(ierr); 2974 } 2975 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2976 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2977 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 2978 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 2979 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 2980 ierr = PetscFree2(perms,flips);CHKERRQ(ierr); 2981 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 2982 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 2983 PetscFunctionReturn(0); 2984 } 2985 2986 /* 2987 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2988 * 2989 * for each coarse dof \phi^c_i: 2990 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2991 * for each fine dof \phi^f_j; 2992 * a_{i,j} = 0; 2993 * for each fine dof \phi^f_k: 2994 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2995 * [^^^ this is = \phi^c_i ^^^] 2996 */ 2997 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2998 { 2999 PetscDS ds; 3000 PetscSection section, cSection; 3001 DMLabel canonical, depth; 3002 Mat cMat, mat; 3003 PetscInt *nnz; 3004 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 3005 PetscInt m, n; 3006 PetscScalar *pointScalar; 3007 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 3008 PetscErrorCode ierr; 3009 3010 PetscFunctionBegin; 3011 ierr = DMGetDefaultSection(refTree,§ion);CHKERRQ(ierr); 3012 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr); 3013 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr); 3014 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr); 3015 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3016 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3017 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr); 3018 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr); 3019 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr); 3020 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr); 3021 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr); 3022 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr); 3023 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */ 3024 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3025 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 3026 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3027 const PetscInt *children; 3028 PetscInt numChildren; 3029 PetscInt i, numChildDof, numSelfDof; 3030 3031 if (canonical) { 3032 PetscInt pCanonical; 3033 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3034 if (p != pCanonical) continue; 3035 } 3036 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3037 if (!numChildren) continue; 3038 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3039 PetscInt child = children[i]; 3040 PetscInt dof; 3041 3042 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3043 numChildDof += dof; 3044 } 3045 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3046 if (!numChildDof || !numSelfDof) continue; 3047 for (f = 0; f < numFields; f++) { 3048 PetscInt selfOff; 3049 3050 if (numSecFields) { /* count the dofs for just this field */ 3051 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3052 PetscInt child = children[i]; 3053 PetscInt dof; 3054 3055 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3056 numChildDof += dof; 3057 } 3058 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3059 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3060 } 3061 else { 3062 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3063 } 3064 for (i = 0; i < numSelfDof; i++) { 3065 nnz[selfOff + i] = numChildDof; 3066 } 3067 } 3068 } 3069 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr); 3070 ierr = PetscFree(nnz);CHKERRQ(ierr); 3071 /* Setp 2: compute entries */ 3072 for (p = pStart; p < pEnd; p++) { 3073 const PetscInt *children; 3074 PetscInt numChildren; 3075 PetscInt i, numChildDof, numSelfDof; 3076 3077 /* same conditions about when entries occur */ 3078 if (canonical) { 3079 PetscInt pCanonical; 3080 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3081 if (p != pCanonical) continue; 3082 } 3083 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3084 if (!numChildren) continue; 3085 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3086 PetscInt child = children[i]; 3087 PetscInt dof; 3088 3089 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3090 numChildDof += dof; 3091 } 3092 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3093 if (!numChildDof || !numSelfDof) continue; 3094 3095 for (f = 0; f < numFields; f++) { 3096 PetscInt selfOff, fComp, numSelfShapes, numChildShapes, parentCell; 3097 PetscInt cellShapeOff; 3098 PetscObject disc; 3099 PetscDualSpace dsp; 3100 PetscClassId classId; 3101 PetscScalar *pointMat; 3102 PetscInt *matRows, *matCols; 3103 PetscInt pO = PETSC_MIN_INT; 3104 const PetscInt *depthNumDof; 3105 3106 if (numSecFields) { 3107 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3108 PetscInt child = children[i]; 3109 PetscInt dof; 3110 3111 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3112 numChildDof += dof; 3113 } 3114 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3115 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3116 ierr = PetscSectionGetFieldComponents(section,f,&fComp);CHKERRQ(ierr); 3117 } 3118 else { 3119 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3120 fComp = 1; 3121 } 3122 numSelfShapes = numSelfDof / fComp; 3123 numChildShapes = numChildDof / fComp; 3124 3125 /* find a cell whose closure contains p */ 3126 if (p >= cStart && p < cEnd) { 3127 parentCell = p; 3128 } 3129 else { 3130 PetscInt *star = NULL; 3131 PetscInt numStar; 3132 3133 parentCell = -1; 3134 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3135 for (i = numStar - 1; i >= 0; i--) { 3136 PetscInt c = star[2 * i]; 3137 3138 if (c >= cStart && c < cEnd) { 3139 parentCell = c; 3140 break; 3141 } 3142 } 3143 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3144 } 3145 /* determine the offset of p's shape functions withing parentCell's shape functions */ 3146 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3147 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr); 3148 if (classId == PETSCFE_CLASSID) { 3149 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr); 3150 } 3151 else if (classId == PETSCFV_CLASSID) { 3152 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr); 3153 } 3154 else { 3155 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");CHKERRQ(ierr); 3156 } 3157 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr); 3158 { 3159 PetscInt *closure = NULL; 3160 PetscInt numClosure; 3161 3162 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3163 for (i = 0, cellShapeOff = 0; i < numClosure; i++) { 3164 PetscInt point = closure[2 * i], pointDepth; 3165 3166 pO = closure[2 * i + 1]; 3167 if (point == p) break; 3168 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3169 cellShapeOff += depthNumDof[pointDepth]; 3170 } 3171 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3172 } 3173 3174 ierr = DMGetWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3175 ierr = DMGetWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3176 matCols = matRows + numSelfShapes; 3177 for (i = 0; i < numSelfShapes; i++) { 3178 matRows[i] = selfOff + i * fComp; 3179 } 3180 { 3181 PetscInt colOff = 0; 3182 3183 for (i = 0; i < numChildren; i++) { 3184 PetscInt child = children[i]; 3185 PetscInt dof, off, j; 3186 3187 if (numSecFields) { 3188 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr); 3189 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr); 3190 } 3191 else { 3192 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr); 3193 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr); 3194 } 3195 3196 for (j = 0; j < dof / fComp; j++) { 3197 matCols[colOff++] = off + j * fComp; 3198 } 3199 } 3200 } 3201 if (classId == PETSCFE_CLASSID) { 3202 PetscFE fe = (PetscFE) disc; 3203 PetscInt fSize; 3204 3205 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 3206 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr); 3207 for (i = 0; i < numSelfShapes; i++) { /* for every shape function */ 3208 PetscQuadrature q; 3209 PetscInt dim, numPoints, j, k; 3210 const PetscReal *points; 3211 const PetscReal *weights; 3212 PetscInt *closure = NULL; 3213 PetscInt numClosure; 3214 PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfShapes - 1 - i) : i); 3215 PetscReal *Bparent; 3216 3217 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr); 3218 ierr = PetscQuadratureGetData(q,&dim,&numPoints,&points,&weights);CHKERRQ(ierr); 3219 ierr = PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3220 for (k = 0; k < numChildShapes; k++) { 3221 pointMat[numChildShapes * i + k] = 0.; 3222 } 3223 for (j = 0; j < numPoints; j++) { 3224 PetscInt childCell = -1; 3225 PetscReal parentValAtPoint; 3226 const PetscReal *pointReal = &points[dim * j]; 3227 const PetscScalar *point; 3228 PetscReal *Bchild; 3229 PetscInt childCellShapeOff, pointMatOff; 3230 #if defined(PETSC_USE_COMPLEX) 3231 PetscInt d; 3232 3233 for (d = 0; d < dim; d++) { 3234 pointScalar[d] = points[dim * j + d]; 3235 } 3236 point = pointScalar; 3237 #else 3238 point = pointReal; 3239 #endif 3240 3241 parentValAtPoint = Bparent[(fSize * j + parentCellShapeDof) * fComp]; 3242 3243 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3244 PetscInt child = children[k]; 3245 PetscInt *star = NULL; 3246 PetscInt numStar, s; 3247 3248 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3249 for (s = numStar - 1; s >= 0; s--) { 3250 PetscInt c = star[2 * s]; 3251 3252 if (c < cStart || c >= cEnd) continue; 3253 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr); 3254 if (childCell >= 0) break; 3255 } 3256 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3257 if (childCell >= 0) break; 3258 } 3259 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");CHKERRQ(ierr); 3260 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 3261 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr); 3262 CoordinatesRefToReal(dim, dim, v0parent, Jparent, pointReal, vtmp); 3263 CoordinatesRealToRef(dim, dim, v0, invJ, vtmp, pointRef); 3264 3265 ierr = PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3266 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3267 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3268 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3269 PetscInt l; 3270 3271 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr); 3272 childDof = depthNumDof[childDepth]; 3273 for (l = 0, childCellShapeOff = 0; l < numClosure; l++) { 3274 PetscInt point = closure[2 * l]; 3275 PetscInt pointDepth; 3276 3277 childO = closure[2 * l + 1]; 3278 if (point == child) break; 3279 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3280 childCellShapeOff += depthNumDof[pointDepth]; 3281 } 3282 if (l == numClosure) { 3283 pointMatOff += childDof; 3284 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3285 } 3286 for (l = 0; l < childDof; l++) { 3287 PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l); 3288 PetscReal childValAtPoint = Bchild[childCellDof * fComp]; 3289 3290 pointMat[i * numChildShapes + pointMatOff + l] += weights[j] * parentValAtPoint * childValAtPoint; 3291 } 3292 pointMatOff += childDof; 3293 } 3294 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3295 ierr = PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3296 } 3297 ierr = PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); 3298 } 3299 } 3300 else { /* just the volume-weighted averages of the children */ 3301 PetscInt childShapeOff; 3302 PetscReal parentVol; 3303 3304 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr); 3305 for (i = 0, childShapeOff = 0; i < numChildren; i++) { 3306 PetscInt child = children[i]; 3307 PetscReal childVol; 3308 3309 if (child < cStart || child >= cEnd) continue; 3310 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr); 3311 pointMat[childShapeOff] = childVol / parentVol; 3312 childShapeOff++; 3313 } 3314 } 3315 /* Insert pointMat into mat */ 3316 for (i = 0; i < fComp; i++) { 3317 PetscInt j; 3318 ierr = MatSetValues(mat,numSelfShapes,matRows,numChildShapes,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr); 3319 3320 for (j = 0; j < numSelfShapes; j++) { 3321 matRows[j]++; 3322 } 3323 for (j = 0; j < numChildShapes; j++) { 3324 matCols[j]++; 3325 } 3326 } 3327 ierr = DMRestoreWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3328 ierr = DMRestoreWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3329 } 3330 } 3331 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr); 3332 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr); 3333 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3334 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3335 *inj = mat; 3336 PetscFunctionReturn(0); 3337 } 3338 3339 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3340 { 3341 PetscDS ds; 3342 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3343 PetscScalar ***refPointFieldMats; 3344 PetscSection refConSec, refSection; 3345 PetscErrorCode ierr; 3346 3347 PetscFunctionBegin; 3348 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3349 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3350 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3351 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3352 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3353 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 3354 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 3355 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 3356 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr); 3357 for (p = pRefStart; p < pRefEnd; p++) { 3358 PetscInt parent, pDof, parentDof; 3359 3360 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3361 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3362 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3363 if (!pDof || !parentDof || parent == p) continue; 3364 3365 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 3366 for (f = 0; f < numFields; f++) { 3367 PetscInt cDof, cOff, numCols, r, fComp; 3368 PetscObject disc; 3369 PetscClassId id; 3370 PetscFE fe = NULL; 3371 PetscFV fv = NULL; 3372 3373 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3374 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 3375 if (id == PETSCFE_CLASSID) { 3376 fe = (PetscFE) disc; 3377 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 3378 } 3379 else if (id == PETSCFV_CLASSID) { 3380 fv = (PetscFV) disc; 3381 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 3382 } 3383 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 3384 3385 if (numFields > 1) { 3386 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3387 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 3388 } 3389 else { 3390 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3391 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 3392 } 3393 3394 for (r = 0; r < cDof; r++) { 3395 rows[r] = cOff + r; 3396 } 3397 numCols = 0; 3398 { 3399 PetscInt aDof, aOff, j; 3400 3401 if (numFields > 1) { 3402 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr); 3403 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr); 3404 } 3405 else { 3406 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr); 3407 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr); 3408 } 3409 3410 for (j = 0; j < aDof; j++) { 3411 cols[numCols++] = aOff + j; 3412 } 3413 } 3414 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3415 /* transpose of constraint matrix */ 3416 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3417 } 3418 } 3419 *childrenMats = refPointFieldMats; 3420 ierr = PetscFree(rows);CHKERRQ(ierr); 3421 ierr = PetscFree(cols);CHKERRQ(ierr); 3422 PetscFunctionReturn(0); 3423 } 3424 3425 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3426 { 3427 PetscDS ds; 3428 PetscScalar ***refPointFieldMats; 3429 PetscInt numFields, pRefStart, pRefEnd, p, f; 3430 PetscSection refConSec, refSection; 3431 PetscErrorCode ierr; 3432 3433 PetscFunctionBegin; 3434 refPointFieldMats = *childrenMats; 3435 *childrenMats = NULL; 3436 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3437 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3438 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3439 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3440 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3441 for (p = pRefStart; p < pRefEnd; p++) { 3442 PetscInt parent, pDof, parentDof; 3443 3444 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3445 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3446 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3447 if (!pDof || !parentDof || parent == p) continue; 3448 3449 for (f = 0; f < numFields; f++) { 3450 PetscInt cDof; 3451 3452 if (numFields > 1) { 3453 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3454 } 3455 else { 3456 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3457 } 3458 3459 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 3460 } 3461 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 3462 } 3463 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 3464 PetscFunctionReturn(0); 3465 } 3466 3467 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3468 { 3469 Mat cMatRef; 3470 PetscObject injRefObj; 3471 PetscErrorCode ierr; 3472 3473 PetscFunctionBegin; 3474 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr); 3475 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr); 3476 *injRef = (Mat) injRefObj; 3477 if (!*injRef) { 3478 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr); 3479 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr); 3480 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3481 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr); 3482 } 3483 PetscFunctionReturn(0); 3484 } 3485 3486 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) 3487 { 3488 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3489 PetscSection globalCoarse, globalFine; 3490 PetscSection localCoarse, localFine, leafIndicesSec; 3491 PetscSection multiRootSec, rootIndicesSec; 3492 PetscInt *leafInds, *rootInds = NULL; 3493 const PetscInt *rootDegrees; 3494 PetscScalar *leafVals = NULL, *rootVals = NULL; 3495 PetscSF coarseToFineEmbedded; 3496 PetscErrorCode ierr; 3497 3498 PetscFunctionBegin; 3499 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3500 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3501 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3502 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3503 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 3504 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr); 3505 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr); 3506 { /* winnow fine points that don't have global dofs out of the sf */ 3507 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3508 const PetscInt *leaves; 3509 3510 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3511 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3512 p = leaves ? leaves[l] : l; 3513 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3514 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3515 if ((dof - cdof) > 0) { 3516 numPointsWithDofs++; 3517 3518 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr); 3519 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr); 3520 } 3521 } 3522 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3523 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr); 3524 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr); 3525 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr); 3526 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);} 3527 for (l = 0, offset = 0; l < nleaves; l++) { 3528 p = leaves ? leaves[l] : l; 3529 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3530 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3531 if ((dof - cdof) > 0) { 3532 PetscInt off, gOff; 3533 PetscInt *pInd; 3534 PetscScalar *pVal = NULL; 3535 3536 pointsWithDofs[offset++] = l; 3537 3538 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3539 3540 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3541 if (gatheredValues) { 3542 PetscInt i; 3543 3544 pVal = &leafVals[off + 1]; 3545 for (i = 0; i < dof; i++) pVal[i] = 0.; 3546 } 3547 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 3548 3549 offsets[0] = 0; 3550 if (numFields) { 3551 PetscInt f; 3552 3553 for (f = 0; f < numFields; f++) { 3554 PetscInt fDof; 3555 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr); 3556 offsets[f + 1] = fDof + offsets[f]; 3557 } 3558 DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 3559 } 3560 else { 3561 DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 3562 } 3563 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);} 3564 } 3565 } 3566 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 3567 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3568 } 3569 3570 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3571 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3572 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3573 3574 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3575 MPI_Datatype threeInt; 3576 PetscMPIInt rank; 3577 PetscInt (*parentNodeAndIdCoarse)[3]; 3578 PetscInt (*parentNodeAndIdFine)[3]; 3579 PetscInt p, nleaves, nleavesToParents; 3580 PetscSF pointSF, sfToParents; 3581 const PetscInt *ilocal; 3582 const PetscSFNode *iremote; 3583 PetscSFNode *iremoteToParents; 3584 PetscInt *ilocalToParents; 3585 3586 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr); 3587 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr); 3588 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr); 3589 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr); 3590 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr); 3591 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 3592 for (p = pStartC; p < pEndC; p++) { 3593 PetscInt parent, childId; 3594 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr); 3595 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3596 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3597 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3598 if (nleaves > 0) { 3599 PetscInt leaf = -1; 3600 3601 if (ilocal) { 3602 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr); 3603 } 3604 else { 3605 leaf = p - pStartC; 3606 } 3607 if (leaf >= 0) { 3608 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3609 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3610 } 3611 } 3612 } 3613 for (p = pStartF; p < pEndF; p++) { 3614 parentNodeAndIdFine[p - pStartF][0] = -1; 3615 parentNodeAndIdFine[p - pStartF][1] = -1; 3616 parentNodeAndIdFine[p - pStartF][2] = -1; 3617 } 3618 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3619 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3620 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3621 PetscInt dof; 3622 3623 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr); 3624 if (dof) { 3625 PetscInt off; 3626 3627 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3628 if (gatheredIndices) { 3629 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3630 } else if (gatheredValues) { 3631 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3632 } 3633 } 3634 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3635 nleavesToParents++; 3636 } 3637 } 3638 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr); 3639 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr); 3640 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3641 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3642 ilocalToParents[nleavesToParents] = p - pStartF; 3643 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3644 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3645 nleavesToParents++; 3646 } 3647 } 3648 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr); 3649 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr); 3650 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3651 3652 coarseToFineEmbedded = sfToParents; 3653 3654 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3655 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr); 3656 } 3657 3658 { /* winnow out coarse points that don't have dofs */ 3659 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3660 PetscSF sfDofsOnly; 3661 3662 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3663 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3664 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3665 if ((dof - cdof) > 0) { 3666 numPointsWithDofs++; 3667 } 3668 } 3669 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3670 for (p = pStartC, offset = 0; p < pEndC; p++) { 3671 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3672 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3673 if ((dof - cdof) > 0) { 3674 pointsWithDofs[offset++] = p - pStartC; 3675 } 3676 } 3677 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr); 3678 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3679 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3680 coarseToFineEmbedded = sfDofsOnly; 3681 } 3682 3683 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3684 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3685 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3686 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr); 3687 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr); 3688 for (p = pStartC; p < pEndC; p++) { 3689 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr); 3690 } 3691 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr); 3692 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr); 3693 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 3694 { /* distribute the leaf section */ 3695 PetscSF multi, multiInv, indicesSF; 3696 PetscInt *remoteOffsets, numRootIndices; 3697 3698 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr); 3699 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr); 3700 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr); 3701 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr); 3702 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 3703 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr); 3704 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 3705 if (gatheredIndices) { 3706 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr); 3707 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3708 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3709 } 3710 if (gatheredValues) { 3711 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr); 3712 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3713 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3714 } 3715 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 3716 } 3717 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3718 ierr = PetscFree(leafInds);CHKERRQ(ierr); 3719 ierr = PetscFree(leafVals);CHKERRQ(ierr); 3720 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3721 *rootMultiSec = multiRootSec; 3722 *multiLeafSec = rootIndicesSec; 3723 if (gatheredIndices) *gatheredIndices = rootInds; 3724 if (gatheredValues) *gatheredValues = rootVals; 3725 PetscFunctionReturn(0); 3726 } 3727 3728 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3729 { 3730 DM refTree; 3731 PetscSection multiRootSec, rootIndicesSec; 3732 PetscSection globalCoarse, globalFine; 3733 PetscSection localCoarse, localFine; 3734 PetscSection cSecRef; 3735 PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd; 3736 Mat injRef; 3737 PetscInt numFields, maxDof; 3738 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3739 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3740 PetscLayout rowMap, colMap; 3741 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3742 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3743 PetscErrorCode ierr; 3744 3745 PetscFunctionBegin; 3746 3747 /* get the templates for the fine-to-coarse injection from the reference tree */ 3748 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 3749 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 3750 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3751 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 3752 3753 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3754 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3755 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3756 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 3757 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3758 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3759 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3760 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 3761 { 3762 PetscInt maxFields = PetscMax(1,numFields) + 1; 3763 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 3764 } 3765 3766 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr); 3767 3768 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr); 3769 3770 /* count indices */ 3771 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 3772 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 3773 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 3774 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 3775 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 3776 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr); 3777 for (p = pStartC; p < pEndC; p++) { 3778 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3779 3780 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3781 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3782 if ((dof - cdof) <= 0) continue; 3783 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3784 3785 rowOffsets[0] = 0; 3786 offsetsCopy[0] = 0; 3787 if (numFields) { 3788 PetscInt f; 3789 3790 for (f = 0; f < numFields; f++) { 3791 PetscInt fDof; 3792 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3793 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3794 } 3795 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3796 } 3797 else { 3798 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3799 rowOffsets[1] = offsetsCopy[0]; 3800 } 3801 3802 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3803 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3804 leafEnd = leafStart + numLeaves; 3805 for (l = leafStart; l < leafEnd; l++) { 3806 PetscInt numIndices, childId, offset; 3807 const PetscInt *childIndices; 3808 3809 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3810 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3811 childId = rootIndices[offset++]; 3812 childIndices = &rootIndices[offset]; 3813 numIndices--; 3814 3815 if (childId == -1) { /* equivalent points: scatter */ 3816 PetscInt i; 3817 3818 for (i = 0; i < numIndices; i++) { 3819 PetscInt colIndex = childIndices[i]; 3820 PetscInt rowIndex = parentIndices[i]; 3821 if (rowIndex < 0) continue; 3822 if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3823 if (colIndex >= colStart && colIndex < colEnd) { 3824 nnzD[rowIndex - rowStart] = 1; 3825 } 3826 else { 3827 nnzO[rowIndex - rowStart] = 1; 3828 } 3829 } 3830 } 3831 else { 3832 PetscInt parentId, f, lim; 3833 3834 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3835 3836 lim = PetscMax(1,numFields); 3837 offsets[0] = 0; 3838 if (numFields) { 3839 PetscInt f; 3840 3841 for (f = 0; f < numFields; f++) { 3842 PetscInt fDof; 3843 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3844 3845 offsets[f + 1] = fDof + offsets[f]; 3846 } 3847 } 3848 else { 3849 PetscInt cDof; 3850 3851 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3852 offsets[1] = cDof; 3853 } 3854 for (f = 0; f < lim; f++) { 3855 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3856 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3857 PetscInt i, numD = 0, numO = 0; 3858 3859 for (i = childStart; i < childEnd; i++) { 3860 PetscInt colIndex = childIndices[i]; 3861 3862 if (colIndex < 0) continue; 3863 if (colIndex >= colStart && colIndex < colEnd) { 3864 numD++; 3865 } 3866 else { 3867 numO++; 3868 } 3869 } 3870 for (i = parentStart; i < parentEnd; i++) { 3871 PetscInt rowIndex = parentIndices[i]; 3872 3873 if (rowIndex < 0) continue; 3874 nnzD[rowIndex - rowStart] += numD; 3875 nnzO[rowIndex - rowStart] += numO; 3876 } 3877 } 3878 } 3879 } 3880 } 3881 /* preallocate */ 3882 ierr = MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);CHKERRQ(ierr); 3883 ierr = PetscFree2(nnzD,nnzO);CHKERRQ(ierr); 3884 /* insert values */ 3885 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3886 for (p = pStartC; p < pEndC; p++) { 3887 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3888 3889 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3890 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3891 if ((dof - cdof) <= 0) continue; 3892 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3893 3894 rowOffsets[0] = 0; 3895 offsetsCopy[0] = 0; 3896 if (numFields) { 3897 PetscInt f; 3898 3899 for (f = 0; f < numFields; f++) { 3900 PetscInt fDof; 3901 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3902 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3903 } 3904 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3905 } 3906 else { 3907 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3908 rowOffsets[1] = offsetsCopy[0]; 3909 } 3910 3911 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3912 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3913 leafEnd = leafStart + numLeaves; 3914 for (l = leafStart; l < leafEnd; l++) { 3915 PetscInt numIndices, childId, offset; 3916 const PetscInt *childIndices; 3917 3918 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3919 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3920 childId = rootIndices[offset++]; 3921 childIndices = &rootIndices[offset]; 3922 numIndices--; 3923 3924 if (childId == -1) { /* equivalent points: scatter */ 3925 PetscInt i; 3926 3927 for (i = 0; i < numIndices; i++) { 3928 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr); 3929 } 3930 } 3931 else { 3932 PetscInt parentId, f, lim; 3933 3934 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3935 3936 lim = PetscMax(1,numFields); 3937 offsets[0] = 0; 3938 if (numFields) { 3939 PetscInt f; 3940 3941 for (f = 0; f < numFields; f++) { 3942 PetscInt fDof; 3943 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3944 3945 offsets[f + 1] = fDof + offsets[f]; 3946 } 3947 } 3948 else { 3949 PetscInt cDof; 3950 3951 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3952 offsets[1] = cDof; 3953 } 3954 for (f = 0; f < lim; f++) { 3955 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3956 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3957 const PetscInt *colIndices = &childIndices[offsets[f]]; 3958 3959 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr); 3960 } 3961 } 3962 } 3963 } 3964 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 3965 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 3966 ierr = PetscFree(parentIndices);CHKERRQ(ierr); 3967 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3968 ierr = PetscFree(rootIndices);CHKERRQ(ierr); 3969 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 3970 3971 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3972 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3973 PetscFunctionReturn(0); 3974 } 3975 3976 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3977 { 3978 PetscDS ds; 3979 PetscSF coarseToFineEmbedded; 3980 PetscSection globalCoarse, globalFine; 3981 PetscSection localCoarse, localFine; 3982 PetscSection aSec, cSec; 3983 PetscSection rootValuesSec; 3984 PetscSection leafValuesSec; 3985 PetscScalar *rootValues, *leafValues; 3986 IS aIS; 3987 const PetscInt *anchors; 3988 Mat cMat; 3989 PetscInt numFields; 3990 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior; 3991 PetscInt aStart, aEnd, cStart, cEnd; 3992 PetscInt *maxChildIds; 3993 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3994 PetscFV fv = NULL; 3995 PetscInt dim, numFVcomps = -1, fvField = -1; 3996 DM cellDM = NULL, gradDM = NULL; 3997 const PetscScalar *cellGeomArray = NULL; 3998 const PetscScalar *gradArray = NULL; 3999 PetscErrorCode ierr; 4000 4001 PetscFunctionBegin; 4002 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4003 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4004 ierr = DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4005 ierr = DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4006 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4007 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4008 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4009 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr); 4010 { /* winnow fine points that don't have global dofs out of the sf */ 4011 PetscInt nleaves, l; 4012 const PetscInt *leaves; 4013 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 4014 4015 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 4016 4017 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 4018 PetscInt p = leaves ? leaves[l] : l; 4019 4020 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4021 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4022 if ((dof - cdof) > 0) { 4023 numPointsWithDofs++; 4024 } 4025 } 4026 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 4027 for (l = 0, offset = 0; l < nleaves; l++) { 4028 PetscInt p = leaves ? leaves[l] : l; 4029 4030 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4031 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4032 if ((dof - cdof) > 0) { 4033 pointsWithDofs[offset++] = l; 4034 } 4035 } 4036 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 4037 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 4038 } 4039 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4040 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 4041 for (p = pStartC; p < pEndC; p++) { 4042 maxChildIds[p - pStartC] = -2; 4043 } 4044 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4045 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4046 4047 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4048 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4049 4050 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 4051 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 4052 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 4053 4054 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 4055 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 4056 4057 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4058 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr); 4059 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr); 4060 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 4061 { 4062 PetscInt maxFields = PetscMax(1,numFields) + 1; 4063 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 4064 } 4065 if (grad) { 4066 PetscInt i; 4067 4068 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr); 4069 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr); 4070 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr); 4071 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr); 4072 for (i = 0; i < PetscMax(1,numFields); i++) { 4073 PetscObject obj; 4074 PetscClassId id; 4075 4076 ierr = DMGetField(coarse, i, &obj);CHKERRQ(ierr); 4077 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr); 4078 if (id == PETSCFV_CLASSID) { 4079 fv = (PetscFV) obj; 4080 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr); 4081 fvField = i; 4082 break; 4083 } 4084 } 4085 } 4086 4087 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4088 PetscInt dof; 4089 PetscInt maxChildId = maxChildIds[p - pStartC]; 4090 PetscInt numValues = 0; 4091 4092 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4093 if (dof < 0) { 4094 dof = -(dof + 1); 4095 } 4096 offsets[0] = 0; 4097 newOffsets[0] = 0; 4098 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4099 PetscInt *closure = NULL, closureSize, cl; 4100 4101 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4102 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4103 PetscInt c = closure[2 * cl], clDof; 4104 4105 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 4106 numValues += clDof; 4107 } 4108 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4109 } 4110 else if (maxChildId == -1) { 4111 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr); 4112 } 4113 /* we will pack the column indices with the field offsets */ 4114 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4115 /* also send the centroid, and the gradient */ 4116 numValues += dim * (1 + numFVcomps); 4117 } 4118 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr); 4119 } 4120 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr); 4121 { 4122 PetscInt numRootValues; 4123 const PetscScalar *coarseArray; 4124 4125 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr); 4126 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr); 4127 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4128 for (p = pStartC; p < pEndC; p++) { 4129 PetscInt numValues; 4130 PetscInt pValOff; 4131 PetscScalar *pVal; 4132 PetscInt maxChildId = maxChildIds[p - pStartC]; 4133 4134 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr); 4135 if (!numValues) { 4136 continue; 4137 } 4138 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr); 4139 pVal = &(rootValues[pValOff]); 4140 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4141 PetscInt closureSize = numValues; 4142 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr); 4143 if (grad && p >= cellStart && p < cellEnd) { 4144 PetscFVCellGeom *cg; 4145 PetscScalar *gradVals = NULL; 4146 PetscInt i; 4147 4148 pVal += (numValues - dim * (1 + numFVcomps)); 4149 4150 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr); 4151 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4152 pVal += dim; 4153 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr); 4154 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4155 } 4156 } 4157 else if (maxChildId == -1) { 4158 PetscInt lDof, lOff, i; 4159 4160 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr); 4161 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr); 4162 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4163 } 4164 } 4165 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4166 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 4167 } 4168 { 4169 PetscSF valuesSF; 4170 PetscInt *remoteOffsetsValues, numLeafValues; 4171 4172 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr); 4173 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr); 4174 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr); 4175 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 4176 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr); 4177 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr); 4178 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr); 4179 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4180 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4181 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr); 4182 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4183 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr); 4184 } 4185 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4186 { 4187 PetscInt maxDof; 4188 PetscInt *rowIndices; 4189 DM refTree; 4190 PetscInt **refPointFieldN; 4191 PetscScalar ***refPointFieldMats; 4192 PetscSection refConSec, refAnSec; 4193 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4194 PetscScalar *pointWork; 4195 4196 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 4197 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4198 ierr = DMGetWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4199 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 4200 ierr = DMGetDS(fine,&ds);CHKERRQ(ierr); 4201 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4202 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4203 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 4204 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 4205 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4206 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 4207 ierr = DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4208 ierr = DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4209 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4210 for (p = leafStart; p < leafEnd; p++) { 4211 PetscInt gDof, gcDof, gOff, lDof; 4212 PetscInt numValues, pValOff; 4213 PetscInt childId; 4214 const PetscScalar *pVal; 4215 const PetscScalar *fvGradData = NULL; 4216 4217 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 4218 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr); 4219 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 4220 if ((gDof - gcDof) <= 0) { 4221 continue; 4222 } 4223 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 4224 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr); 4225 if (!numValues) continue; 4226 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr); 4227 pVal = &leafValues[pValOff]; 4228 offsets[0] = 0; 4229 offsetsCopy[0] = 0; 4230 newOffsets[0] = 0; 4231 newOffsetsCopy[0] = 0; 4232 childId = cids[p - pStartF]; 4233 if (numFields) { 4234 PetscInt f; 4235 for (f = 0; f < numFields; f++) { 4236 PetscInt rowDof; 4237 4238 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 4239 offsets[f + 1] = offsets[f] + rowDof; 4240 offsetsCopy[f + 1] = offsets[f + 1]; 4241 /* TODO: closure indices */ 4242 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4243 } 4244 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 4245 } 4246 else { 4247 offsets[0] = 0; 4248 offsets[1] = lDof; 4249 newOffsets[0] = 0; 4250 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4251 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 4252 } 4253 if (childId == -1) { /* no child interpolation: one nnz per */ 4254 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr); 4255 } else { 4256 PetscInt f; 4257 4258 if (grad && p >= cellStart && p < cellEnd) { 4259 numValues -= (dim * (1 + numFVcomps)); 4260 fvGradData = &pVal[numValues]; 4261 } 4262 for (f = 0; f < PetscMax(1,numFields); f++) { 4263 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4264 PetscInt numRows = offsets[f+1] - offsets[f]; 4265 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4266 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4267 PetscScalar *rVal = &pointWork[offsets[f]]; 4268 PetscInt i, j; 4269 4270 #if 0 4271 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); 4272 #endif 4273 for (i = 0; i < numRows; i++) { 4274 PetscScalar val = 0.; 4275 for (j = 0; j < numCols; j++) { 4276 val += childMat[i * numCols + j] * cVal[j]; 4277 } 4278 rVal[i] = val; 4279 } 4280 if (f == fvField && p >= cellStart && p < cellEnd) { 4281 PetscReal centroid[3]; 4282 PetscScalar diff[3]; 4283 const PetscScalar *parentCentroid = &fvGradData[0]; 4284 const PetscScalar *gradient = &fvGradData[dim]; 4285 4286 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr); 4287 for (i = 0; i < dim; i++) { 4288 diff[i] = centroid[i] - parentCentroid[i]; 4289 } 4290 for (i = 0; i < numFVcomps; i++) { 4291 PetscScalar val = 0.; 4292 4293 for (j = 0; j < dim; j++) { 4294 val += gradient[dim * i + j] * diff[j]; 4295 } 4296 rVal[i] += val; 4297 } 4298 } 4299 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr); 4300 } 4301 } 4302 } 4303 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4304 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4305 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4306 } 4307 ierr = PetscFree(leafValues);CHKERRQ(ierr); 4308 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr); 4309 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 4310 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 4311 PetscFunctionReturn(0); 4312 } 4313 4314 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4315 { 4316 PetscDS ds; 4317 DM refTree; 4318 PetscSection multiRootSec, rootIndicesSec; 4319 PetscSection globalCoarse, globalFine; 4320 PetscSection localCoarse, localFine; 4321 PetscSection cSecRef; 4322 PetscInt *parentIndices, pRefStart, pRefEnd; 4323 PetscScalar *rootValues, *parentValues; 4324 Mat injRef; 4325 PetscInt numFields, maxDof; 4326 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4327 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4328 PetscLayout rowMap, colMap; 4329 PetscInt rowStart, rowEnd, colStart, colEnd; 4330 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4331 PetscErrorCode ierr; 4332 4333 PetscFunctionBegin; 4334 4335 /* get the templates for the fine-to-coarse injection from the reference tree */ 4336 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4337 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4338 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 4339 ierr = DMGetDS(coarse,&ds);CHKERRQ(ierr); 4340 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4341 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 4342 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4343 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 4344 4345 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4346 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4347 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4348 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 4349 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4350 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4351 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4352 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 4353 { 4354 PetscInt maxFields = PetscMax(1,numFields) + 1; 4355 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 4356 } 4357 4358 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr); 4359 4360 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr); 4361 4362 /* count indices */ 4363 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr); 4364 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr); 4365 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 4366 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 4367 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 4368 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 4369 /* insert values */ 4370 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4371 for (p = pStartC; p < pEndC; p++) { 4372 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4373 PetscBool contribute = PETSC_FALSE; 4374 4375 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4376 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 4377 if ((dof - cdof) <= 0) continue; 4378 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr); 4379 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 4380 4381 rowOffsets[0] = 0; 4382 offsetsCopy[0] = 0; 4383 if (numFields) { 4384 PetscInt f; 4385 4386 for (f = 0; f < numFields; f++) { 4387 PetscInt fDof; 4388 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 4389 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4390 } 4391 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 4392 } 4393 else { 4394 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 4395 rowOffsets[1] = offsetsCopy[0]; 4396 } 4397 4398 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 4399 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 4400 leafEnd = leafStart + numLeaves; 4401 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4402 for (l = leafStart; l < leafEnd; l++) { 4403 PetscInt numIndices, childId, offset; 4404 const PetscScalar *childValues; 4405 4406 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 4407 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 4408 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4409 childValues = &rootValues[offset]; 4410 numIndices--; 4411 4412 if (childId == -2) { /* skip */ 4413 continue; 4414 } else if (childId == -1) { /* equivalent points: scatter */ 4415 PetscInt m; 4416 4417 contribute = PETSC_TRUE; 4418 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4419 } else { /* contributions from children: sum with injectors from reference tree */ 4420 PetscInt parentId, f, lim; 4421 4422 contribute = PETSC_TRUE; 4423 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 4424 4425 lim = PetscMax(1,numFields); 4426 offsets[0] = 0; 4427 if (numFields) { 4428 PetscInt f; 4429 4430 for (f = 0; f < numFields; f++) { 4431 PetscInt fDof; 4432 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4433 4434 offsets[f + 1] = fDof + offsets[f]; 4435 } 4436 } 4437 else { 4438 PetscInt cDof; 4439 4440 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4441 offsets[1] = cDof; 4442 } 4443 for (f = 0; f < lim; f++) { 4444 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4445 PetscInt n = offsets[f+1]-offsets[f]; 4446 PetscInt m = rowOffsets[f+1]-rowOffsets[f]; 4447 PetscInt i, j; 4448 const PetscScalar *colValues = &childValues[offsets[f]]; 4449 4450 for (i = 0; i < m; i++) { 4451 PetscScalar val = 0.; 4452 for (j = 0; j < n; j++) { 4453 val += childMat[n * i + j] * colValues[j]; 4454 } 4455 parentValues[rowOffsets[f] + i] += val; 4456 } 4457 } 4458 } 4459 } 4460 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);} 4461 } 4462 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4463 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4464 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr); 4465 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4466 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4467 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4468 PetscFunctionReturn(0); 4469 } 4470 4471 /*@ 4472 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4473 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4474 coarsening and refinement at the same time. 4475 4476 collective 4477 4478 Input Parameters: 4479 + dmIn - The DMPlex mesh for the input vector 4480 . vecIn - The input vector 4481 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4482 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4483 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4484 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4485 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4486 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4487 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4488 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4489 point j in dmOut is not a leaf of sfRefine. 4490 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4491 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4492 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4493 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4494 - time - Used if boundary values are time dependent. 4495 4496 Output Parameters: 4497 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered 4498 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4499 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4500 coarse points to fine points. 4501 4502 Level: developer 4503 4504 .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients() 4505 @*/ 4506 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4507 { 4508 PetscErrorCode ierr; 4509 4510 PetscFunctionBegin; 4511 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4512 if (sfRefine) { 4513 Vec vecInLocal; 4514 DM dmGrad = NULL; 4515 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4516 4517 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4518 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4519 { 4520 PetscInt numFields, i; 4521 4522 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4523 for (i = 0; i < numFields; i++) { 4524 PetscObject obj; 4525 PetscClassId classid; 4526 4527 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4528 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4529 if (classid == PETSCFV_CLASSID) { 4530 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4531 break; 4532 } 4533 } 4534 } 4535 if (useBCs) { 4536 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4537 } 4538 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4539 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4540 if (dmGrad) { 4541 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4542 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4543 } 4544 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr); 4545 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4546 if (dmGrad) { 4547 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4548 } 4549 } 4550 if (sfCoarsen) { 4551 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr); 4552 } 4553 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4554 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4555 PetscFunctionReturn(0); 4556 } 4557