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