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