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