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