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 = PetscCalloc2(maxFields,&perms,maxFields,&flips);CHKERRQ(ierr); 2282 2283 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2284 PetscInt dof, matSize = 0; 2285 PetscInt aDof = 0; 2286 PetscInt cDof = 0; 2287 PetscInt maxChildId = maxChildIds[p - pStartC]; 2288 PetscInt numRowIndices = 0; 2289 PetscInt numColIndices = 0; 2290 PetscInt f; 2291 2292 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2293 if (dof < 0) { 2294 dof = -(dof + 1); 2295 } 2296 if (p >= aStart && p < aEnd) { 2297 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2298 } 2299 if (p >= cStart && p < cEnd) { 2300 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr); 2301 } 2302 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2303 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2304 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2305 PetscInt *closure = NULL, closureSize, cl; 2306 2307 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2308 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2309 PetscInt c = closure[2 * cl], clDof; 2310 2311 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 2312 numRowIndices += clDof; 2313 for (f = 0; f < numFields; f++) { 2314 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr); 2315 offsets[f + 1] += clDof; 2316 } 2317 } 2318 for (f = 0; f < numFields; f++) { 2319 offsets[f + 1] += offsets[f]; 2320 newOffsets[f + 1] = offsets[f + 1]; 2321 } 2322 /* get the number of indices needed and their field offsets */ 2323 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2324 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2325 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2326 numColIndices = numRowIndices; 2327 matSize = 0; 2328 } 2329 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2330 matSize = 0; 2331 for (f = 0; f < numFields; f++) { 2332 PetscInt numRow, numCol; 2333 2334 numRow = offsets[f + 1] - offsets[f]; 2335 numCol = newOffsets[f + 1] - newOffsets[f]; 2336 matSize += numRow * numCol; 2337 } 2338 } 2339 else { 2340 matSize = numRowIndices * numColIndices; 2341 } 2342 } else if (maxChildId == -1) { 2343 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2344 PetscInt aOff, a; 2345 2346 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2347 for (f = 0; f < numFields; f++) { 2348 PetscInt fDof; 2349 2350 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2351 offsets[f+1] = fDof; 2352 } 2353 for (a = 0; a < aDof; a++) { 2354 PetscInt anchor = anchors[a + aOff], aLocalDof; 2355 2356 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr); 2357 numColIndices += aLocalDof; 2358 for (f = 0; f < numFields; f++) { 2359 PetscInt fDof; 2360 2361 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2362 newOffsets[f+1] += fDof; 2363 } 2364 } 2365 if (numFields) { 2366 matSize = 0; 2367 for (f = 0; f < numFields; f++) { 2368 matSize += offsets[f+1] * newOffsets[f+1]; 2369 } 2370 } 2371 else { 2372 matSize = numColIndices * dof; 2373 } 2374 } 2375 else { /* no children, and no constraints on dofs: just get the global indices */ 2376 numColIndices = dof; 2377 matSize = 0; 2378 } 2379 } 2380 /* we will pack the column indices with the field offsets */ 2381 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr); 2382 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr); 2383 } 2384 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr); 2385 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr); 2386 { 2387 PetscInt numRootIndices, numRootMatrices; 2388 2389 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 2390 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr); 2391 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr); 2392 for (p = pStartC; p < pEndC; p++) { 2393 PetscInt numRowIndices, numColIndices, matSize, dof; 2394 PetscInt pIndOff, pMatOff, f; 2395 PetscInt *pInd; 2396 PetscInt maxChildId = maxChildIds[p - pStartC]; 2397 PetscScalar *pMat = NULL; 2398 2399 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2400 if (!numColIndices) { 2401 continue; 2402 } 2403 for (f = 0; f <= numFields; f++) { 2404 offsets[f] = 0; 2405 newOffsets[f] = 0; 2406 offsetsCopy[f] = 0; 2407 newOffsetsCopy[f] = 0; 2408 } 2409 numColIndices -= 2 * numFields; 2410 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2411 pInd = &(rootIndices[pIndOff]); 2412 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr); 2413 if (matSize) { 2414 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2415 pMat = &rootMatrices[pMatOff]; 2416 } 2417 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2418 if (dof < 0) { 2419 dof = -(dof + 1); 2420 } 2421 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2422 PetscInt i, j; 2423 PetscInt numRowIndices = matSize / numColIndices; 2424 2425 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2426 PetscInt numIndices, *indices; 2427 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2428 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2429 for (i = 0; i < numColIndices; i++) { 2430 pInd[i] = indices[i]; 2431 } 2432 for (i = 0; i < numFields; i++) { 2433 pInd[numColIndices + i] = offsets[i+1]; 2434 pInd[numColIndices + numFields + i] = offsets[i+1]; 2435 } 2436 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2437 } 2438 else { 2439 PetscInt closureSize, *closure = NULL, cl; 2440 PetscScalar *pMatIn, *pMatModified; 2441 PetscInt numPoints,*points; 2442 2443 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2444 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2445 for (j = 0; j < numRowIndices; j++) { 2446 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2447 } 2448 } 2449 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2450 for (f = 0; f < maxFields; f++) { 2451 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2452 else {ierr = PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2453 } 2454 if (numFields) { 2455 for (cl = 0; cl < closureSize; cl++) { 2456 PetscInt c = closure[2 * cl]; 2457 2458 for (f = 0; f < numFields; f++) { 2459 PetscInt fDof; 2460 2461 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr); 2462 offsets[f + 1] += fDof; 2463 } 2464 } 2465 for (f = 0; f < numFields; f++) { 2466 offsets[f + 1] += offsets[f]; 2467 newOffsets[f + 1] = offsets[f + 1]; 2468 } 2469 } 2470 /* TODO : flips here ? */ 2471 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2472 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2473 for (f = 0; f < maxFields; f++) { 2474 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2475 else {ierr = PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2476 } 2477 for (f = 0; f < maxFields; f++) { 2478 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2479 else {ierr = PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2480 } 2481 if (!numFields) { 2482 for (i = 0; i < numRowIndices * numColIndices; i++) { 2483 pMat[i] = pMatModified[i]; 2484 } 2485 } 2486 else { 2487 PetscInt i, j, count; 2488 for (f = 0, count = 0; f < numFields; f++) { 2489 for (i = offsets[f]; i < offsets[f+1]; i++) { 2490 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2491 pMat[count] = pMatModified[i * numColIndices + j]; 2492 } 2493 } 2494 } 2495 } 2496 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);CHKERRQ(ierr); 2497 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2498 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2499 if (numFields) { 2500 for (f = 0; f < numFields; f++) { 2501 pInd[numColIndices + f] = offsets[f+1]; 2502 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2503 } 2504 for (cl = 0; cl < numPoints; cl++) { 2505 PetscInt globalOff, c = points[2*cl]; 2506 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2507 DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd); 2508 } 2509 } else { 2510 for (cl = 0; cl < numPoints; cl++) { 2511 PetscInt c = points[2*cl], globalOff; 2512 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2513 2514 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2515 DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd); 2516 } 2517 } 2518 for (f = 0; f < maxFields; f++) { 2519 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2520 else {ierr = PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2521 } 2522 ierr = DMRestoreWorkArray(coarse,numPoints,PETSC_SCALAR,&points);CHKERRQ(ierr); 2523 } 2524 } 2525 else if (matSize) { 2526 PetscInt cOff; 2527 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2528 2529 numRowIndices = matSize / numColIndices; 2530 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2531 ierr = DMGetWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2532 ierr = DMGetWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2533 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr); 2534 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2535 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2536 if (numFields) { 2537 for (f = 0; f < numFields; f++) { 2538 PetscInt fDof; 2539 2540 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr); 2541 offsets[f + 1] = fDof; 2542 for (a = 0; a < aDof; a++) { 2543 PetscInt anchor = anchors[a + aOff]; 2544 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2545 newOffsets[f + 1] += fDof; 2546 } 2547 } 2548 for (f = 0; f < numFields; f++) { 2549 offsets[f + 1] += offsets[f]; 2550 offsetsCopy[f + 1] = offsets[f + 1]; 2551 newOffsets[f + 1] += newOffsets[f]; 2552 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2553 } 2554 DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);CHKERRQ(ierr); 2555 for (a = 0; a < aDof; a++) { 2556 PetscInt anchor = anchors[a + aOff], lOff; 2557 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2558 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);CHKERRQ(ierr); 2559 } 2560 } 2561 else { 2562 DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);CHKERRQ(ierr); 2563 for (a = 0; a < aDof; a++) { 2564 PetscInt anchor = anchors[a + aOff], lOff; 2565 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2566 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);CHKERRQ(ierr); 2567 } 2568 } 2569 if (numFields) { 2570 PetscInt count, a; 2571 2572 for (f = 0, count = 0; f < numFields; f++) { 2573 PetscInt iSize = offsets[f + 1] - offsets[f]; 2574 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2575 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2576 count += iSize * jSize; 2577 pInd[numColIndices + f] = offsets[f+1]; 2578 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2579 } 2580 for (a = 0; a < aDof; a++) { 2581 PetscInt anchor = anchors[a + aOff]; 2582 PetscInt gOff; 2583 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2584 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2585 } 2586 } 2587 else { 2588 PetscInt a; 2589 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2590 for (a = 0; a < aDof; a++) { 2591 PetscInt anchor = anchors[a + aOff]; 2592 PetscInt gOff; 2593 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2594 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2595 } 2596 } 2597 ierr = DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2598 ierr = DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2599 } 2600 else { 2601 PetscInt gOff; 2602 2603 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 2604 if (numFields) { 2605 for (f = 0; f < numFields; f++) { 2606 PetscInt fDof; 2607 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2608 offsets[f + 1] = fDof + offsets[f]; 2609 } 2610 for (f = 0; f < numFields; f++) { 2611 pInd[numColIndices + f] = offsets[f+1]; 2612 pInd[numColIndices + numFields + f] = offsets[f+1]; 2613 } 2614 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2615 } 2616 else { 2617 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2618 } 2619 } 2620 } 2621 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 2622 } 2623 { 2624 PetscSF indicesSF, matricesSF; 2625 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2626 2627 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 2628 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr); 2629 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr); 2630 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr); 2631 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr); 2632 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr); 2633 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 2634 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr); 2635 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr); 2636 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr); 2637 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr); 2638 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr); 2639 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2640 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2641 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2642 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2643 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr); 2644 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 2645 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr); 2646 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 2647 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr); 2648 } 2649 /* count to preallocate */ 2650 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 2651 { 2652 PetscInt nGlobal; 2653 PetscInt *dnnz, *onnz; 2654 PetscLayout rowMap, colMap; 2655 PetscInt rowStart, rowEnd, colStart, colEnd; 2656 PetscInt maxDof; 2657 PetscInt *rowIndices; 2658 DM refTree; 2659 PetscInt **refPointFieldN; 2660 PetscScalar ***refPointFieldMats; 2661 PetscSection refConSec, refAnSec; 2662 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2663 PetscScalar *pointWork; 2664 2665 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr); 2666 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr); 2667 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 2668 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 2669 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 2670 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 2671 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 2672 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 2673 ierr = PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 2674 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2675 for (p = leafStart; p < leafEnd; p++) { 2676 PetscInt gDof, gcDof, gOff; 2677 PetscInt numColIndices, pIndOff, *pInd; 2678 PetscInt matSize; 2679 PetscInt i; 2680 2681 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2682 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2683 if ((gDof - gcDof) <= 0) { 2684 continue; 2685 } 2686 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2687 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2688 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2689 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2690 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2691 numColIndices -= 2 * numFields; 2692 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2693 pInd = &leafIndices[pIndOff]; 2694 offsets[0] = 0; 2695 offsetsCopy[0] = 0; 2696 newOffsets[0] = 0; 2697 newOffsetsCopy[0] = 0; 2698 if (numFields) { 2699 PetscInt f; 2700 for (f = 0; f < numFields; f++) { 2701 PetscInt rowDof; 2702 2703 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2704 offsets[f + 1] = offsets[f] + rowDof; 2705 offsetsCopy[f + 1] = offsets[f + 1]; 2706 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2707 numD[f] = 0; 2708 numO[f] = 0; 2709 } 2710 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2711 for (f = 0; f < numFields; f++) { 2712 PetscInt colOffset = newOffsets[f]; 2713 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2714 2715 for (i = 0; i < numFieldCols; i++) { 2716 PetscInt gInd = pInd[i + colOffset]; 2717 2718 if (gInd >= colStart && gInd < colEnd) { 2719 numD[f]++; 2720 } 2721 else if (gInd >= 0) { /* negative means non-entry */ 2722 numO[f]++; 2723 } 2724 } 2725 } 2726 } 2727 else { 2728 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2729 numD[0] = 0; 2730 numO[0] = 0; 2731 for (i = 0; i < numColIndices; i++) { 2732 PetscInt gInd = pInd[i]; 2733 2734 if (gInd >= colStart && gInd < colEnd) { 2735 numD[0]++; 2736 } 2737 else if (gInd >= 0) { /* negative means non-entry */ 2738 numO[0]++; 2739 } 2740 } 2741 } 2742 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2743 if (!matSize) { /* incoming matrix is identity */ 2744 PetscInt childId; 2745 2746 childId = childIds[p-pStartF]; 2747 if (childId < 0) { /* no child interpolation: one nnz per */ 2748 if (numFields) { 2749 PetscInt f; 2750 for (f = 0; f < numFields; f++) { 2751 PetscInt numRows = offsets[f+1] - offsets[f], row; 2752 for (row = 0; row < numRows; row++) { 2753 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2754 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2755 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2756 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2757 dnnz[gIndFine - rowStart] = 1; 2758 } 2759 else if (gIndCoarse >= 0) { /* remote */ 2760 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2761 onnz[gIndFine - rowStart] = 1; 2762 } 2763 else { /* constrained */ 2764 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2765 } 2766 } 2767 } 2768 } 2769 else { 2770 PetscInt i; 2771 for (i = 0; i < gDof; i++) { 2772 PetscInt gIndCoarse = pInd[i]; 2773 PetscInt gIndFine = rowIndices[i]; 2774 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2775 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2776 dnnz[gIndFine - rowStart] = 1; 2777 } 2778 else if (gIndCoarse >= 0) { /* remote */ 2779 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2780 onnz[gIndFine - rowStart] = 1; 2781 } 2782 else { /* constrained */ 2783 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2784 } 2785 } 2786 } 2787 } 2788 else { /* interpolate from all */ 2789 if (numFields) { 2790 PetscInt f; 2791 for (f = 0; f < numFields; f++) { 2792 PetscInt numRows = offsets[f+1] - offsets[f], row; 2793 for (row = 0; row < numRows; row++) { 2794 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2795 if (gIndFine >= 0) { 2796 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2797 dnnz[gIndFine - rowStart] = numD[f]; 2798 onnz[gIndFine - rowStart] = numO[f]; 2799 } 2800 } 2801 } 2802 } 2803 else { 2804 PetscInt i; 2805 for (i = 0; i < gDof; i++) { 2806 PetscInt gIndFine = rowIndices[i]; 2807 if (gIndFine >= 0) { 2808 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2809 dnnz[gIndFine - rowStart] = numD[0]; 2810 onnz[gIndFine - rowStart] = numO[0]; 2811 } 2812 } 2813 } 2814 } 2815 } 2816 else { /* interpolate from all */ 2817 if (numFields) { 2818 PetscInt f; 2819 for (f = 0; f < numFields; f++) { 2820 PetscInt numRows = offsets[f+1] - offsets[f], row; 2821 for (row = 0; row < numRows; row++) { 2822 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2823 if (gIndFine >= 0) { 2824 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2825 dnnz[gIndFine - rowStart] = numD[f]; 2826 onnz[gIndFine - rowStart] = numO[f]; 2827 } 2828 } 2829 } 2830 } 2831 else { /* every dof get a full row */ 2832 PetscInt i; 2833 for (i = 0; i < gDof; i++) { 2834 PetscInt gIndFine = rowIndices[i]; 2835 if (gIndFine >= 0) { 2836 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2837 dnnz[gIndFine - rowStart] = numD[0]; 2838 onnz[gIndFine - rowStart] = numO[0]; 2839 } 2840 } 2841 } 2842 } 2843 } 2844 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr); 2845 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2846 2847 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 2848 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2849 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 2850 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 2851 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 2852 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr); 2853 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr); 2854 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr); 2855 for (p = leafStart; p < leafEnd; p++) { 2856 PetscInt gDof, gcDof, gOff; 2857 PetscInt numColIndices, pIndOff, *pInd; 2858 PetscInt matSize; 2859 PetscInt childId; 2860 2861 2862 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2863 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2864 if ((gDof - gcDof) <= 0) { 2865 continue; 2866 } 2867 childId = childIds[p-pStartF]; 2868 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2869 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2870 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2871 numColIndices -= 2 * numFields; 2872 pInd = &leafIndices[pIndOff]; 2873 offsets[0] = 0; 2874 offsetsCopy[0] = 0; 2875 newOffsets[0] = 0; 2876 newOffsetsCopy[0] = 0; 2877 rowOffsets[0] = 0; 2878 if (numFields) { 2879 PetscInt f; 2880 for (f = 0; f < numFields; f++) { 2881 PetscInt rowDof; 2882 2883 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2884 offsets[f + 1] = offsets[f] + rowDof; 2885 offsetsCopy[f + 1] = offsets[f + 1]; 2886 rowOffsets[f + 1] = pInd[numColIndices + f]; 2887 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2888 } 2889 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2890 } 2891 else { 2892 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2893 } 2894 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2895 if (!matSize) { /* incoming matrix is identity */ 2896 if (childId < 0) { /* no child interpolation: scatter */ 2897 if (numFields) { 2898 PetscInt f; 2899 for (f = 0; f < numFields; f++) { 2900 PetscInt numRows = offsets[f+1] - offsets[f], row; 2901 for (row = 0; row < numRows; row++) { 2902 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2903 } 2904 } 2905 } 2906 else { 2907 PetscInt numRows = gDof, row; 2908 for (row = 0; row < numRows; row++) { 2909 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2910 } 2911 } 2912 } 2913 else { /* interpolate from all */ 2914 if (numFields) { 2915 PetscInt f; 2916 for (f = 0; f < numFields; f++) { 2917 PetscInt numRows = offsets[f+1] - offsets[f]; 2918 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2919 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2920 } 2921 } 2922 else { 2923 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2924 } 2925 } 2926 } 2927 else { /* interpolate from all */ 2928 PetscInt pMatOff; 2929 PetscScalar *pMat; 2930 2931 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2932 pMat = &leafMatrices[pMatOff]; 2933 if (childId < 0) { /* copy the incoming matrix */ 2934 if (numFields) { 2935 PetscInt f, count; 2936 for (f = 0, count = 0; f < numFields; f++) { 2937 PetscInt numRows = offsets[f+1]-offsets[f]; 2938 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2939 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2940 PetscScalar *inMat = &pMat[count]; 2941 2942 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2943 count += numCols * numInRows; 2944 } 2945 } 2946 else { 2947 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2948 } 2949 } 2950 else { /* multiply the incoming matrix by the child interpolation */ 2951 if (numFields) { 2952 PetscInt f, count; 2953 for (f = 0, count = 0; f < numFields; f++) { 2954 PetscInt numRows = offsets[f+1]-offsets[f]; 2955 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2956 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2957 PetscScalar *inMat = &pMat[count]; 2958 PetscInt i, j, k; 2959 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2960 for (i = 0; i < numRows; i++) { 2961 for (j = 0; j < numCols; j++) { 2962 PetscScalar val = 0.; 2963 for (k = 0; k < numInRows; k++) { 2964 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2965 } 2966 pointWork[i * numCols + j] = val; 2967 } 2968 } 2969 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2970 count += numCols * numInRows; 2971 } 2972 } 2973 else { /* every dof gets a full row */ 2974 PetscInt numRows = gDof; 2975 PetscInt numCols = numColIndices; 2976 PetscInt numInRows = matSize / numColIndices; 2977 PetscInt i, j, k; 2978 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2979 for (i = 0; i < numRows; i++) { 2980 for (j = 0; j < numCols; j++) { 2981 PetscScalar val = 0.; 2982 for (k = 0; k < numInRows; k++) { 2983 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2984 } 2985 pointWork[i * numCols + j] = val; 2986 } 2987 } 2988 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 2989 } 2990 } 2991 } 2992 } 2993 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2994 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2995 ierr = PetscFree(pointWork);CHKERRQ(ierr); 2996 } 2997 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2998 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2999 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3000 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 3001 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 3002 ierr = PetscFree2(perms,flips);CHKERRQ(ierr); 3003 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 3004 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 3005 PetscFunctionReturn(0); 3006 } 3007 3008 #undef __FUNCT__ 3009 #define __FUNCT__ "DMPlexComputeInjectorReferenceTree" 3010 /* 3011 * Assuming a nodal basis (w.r.t. the dual basis) basis: 3012 * 3013 * for each coarse dof \phi^c_i: 3014 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 3015 * for each fine dof \phi^f_j; 3016 * a_{i,j} = 0; 3017 * for each fine dof \phi^f_k: 3018 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 3019 * [^^^ this is = \phi^c_i ^^^] 3020 */ 3021 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 3022 { 3023 PetscDS ds; 3024 PetscSection section, cSection; 3025 DMLabel canonical, depth; 3026 Mat cMat, mat; 3027 PetscInt *nnz; 3028 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 3029 PetscInt m, n; 3030 PetscScalar *pointScalar; 3031 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 3032 PetscErrorCode ierr; 3033 3034 PetscFunctionBegin; 3035 ierr = DMGetDefaultSection(refTree,§ion);CHKERRQ(ierr); 3036 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr); 3037 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr); 3038 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr); 3039 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3040 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3041 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr); 3042 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr); 3043 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr); 3044 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr); 3045 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr); 3046 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr); 3047 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */ 3048 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3049 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 3050 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3051 const PetscInt *children; 3052 PetscInt numChildren; 3053 PetscInt i, numChildDof, numSelfDof; 3054 3055 if (canonical) { 3056 PetscInt pCanonical; 3057 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3058 if (p != pCanonical) continue; 3059 } 3060 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3061 if (!numChildren) continue; 3062 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3063 PetscInt child = children[i]; 3064 PetscInt dof; 3065 3066 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3067 numChildDof += dof; 3068 } 3069 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3070 if (!numChildDof || !numSelfDof) continue; 3071 for (f = 0; f < numFields; f++) { 3072 PetscInt selfOff; 3073 3074 if (numSecFields) { /* count the dofs for just this field */ 3075 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3076 PetscInt child = children[i]; 3077 PetscInt dof; 3078 3079 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3080 numChildDof += dof; 3081 } 3082 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3083 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3084 } 3085 else { 3086 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3087 } 3088 for (i = 0; i < numSelfDof; i++) { 3089 nnz[selfOff + i] = numChildDof; 3090 } 3091 } 3092 } 3093 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr); 3094 ierr = PetscFree(nnz);CHKERRQ(ierr); 3095 /* Setp 2: compute entries */ 3096 for (p = pStart; p < pEnd; p++) { 3097 const PetscInt *children; 3098 PetscInt numChildren; 3099 PetscInt i, numChildDof, numSelfDof; 3100 3101 /* same conditions about when entries occur */ 3102 if (canonical) { 3103 PetscInt pCanonical; 3104 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3105 if (p != pCanonical) continue; 3106 } 3107 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3108 if (!numChildren) continue; 3109 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3110 PetscInt child = children[i]; 3111 PetscInt dof; 3112 3113 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3114 numChildDof += dof; 3115 } 3116 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3117 if (!numChildDof || !numSelfDof) continue; 3118 3119 for (f = 0; f < numFields; f++) { 3120 PetscInt selfOff, fComp, numSelfShapes, numChildShapes, parentCell; 3121 PetscInt cellShapeOff; 3122 PetscObject disc; 3123 PetscDualSpace dsp; 3124 PetscClassId classId; 3125 PetscScalar *pointMat; 3126 PetscInt *matRows, *matCols; 3127 PetscInt pO = PETSC_MIN_INT; 3128 const PetscInt *depthNumDof; 3129 3130 if (numSecFields) { 3131 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3132 PetscInt child = children[i]; 3133 PetscInt dof; 3134 3135 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3136 numChildDof += dof; 3137 } 3138 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3139 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3140 ierr = PetscSectionGetFieldComponents(section,f,&fComp);CHKERRQ(ierr); 3141 } 3142 else { 3143 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3144 fComp = 1; 3145 } 3146 numSelfShapes = numSelfDof / fComp; 3147 numChildShapes = numChildDof / fComp; 3148 3149 /* find a cell whose closure contains p */ 3150 if (p >= cStart && p < cEnd) { 3151 parentCell = p; 3152 } 3153 else { 3154 PetscInt *star = NULL; 3155 PetscInt numStar; 3156 3157 parentCell = -1; 3158 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3159 for (i = numStar - 1; i >= 0; i--) { 3160 PetscInt c = star[2 * i]; 3161 3162 if (c >= cStart && c < cEnd) { 3163 parentCell = c; 3164 break; 3165 } 3166 } 3167 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3168 } 3169 /* determine the offset of p's shape functions withing parentCell's shape functions */ 3170 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3171 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr); 3172 if (classId == PETSCFE_CLASSID) { 3173 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr); 3174 } 3175 else if (classId == PETSCFV_CLASSID) { 3176 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr); 3177 } 3178 else { 3179 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");CHKERRQ(ierr); 3180 } 3181 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr); 3182 { 3183 PetscInt *closure = NULL; 3184 PetscInt numClosure; 3185 3186 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3187 for (i = 0, cellShapeOff = 0; i < numClosure; i++) { 3188 PetscInt point = closure[2 * i], pointDepth; 3189 3190 pO = closure[2 * i + 1]; 3191 if (point == p) break; 3192 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3193 cellShapeOff += depthNumDof[pointDepth]; 3194 } 3195 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3196 } 3197 3198 ierr = DMGetWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3199 ierr = DMGetWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3200 matCols = matRows + numSelfShapes; 3201 for (i = 0; i < numSelfShapes; i++) { 3202 matRows[i] = selfOff + i * fComp; 3203 } 3204 { 3205 PetscInt colOff = 0; 3206 3207 for (i = 0; i < numChildren; i++) { 3208 PetscInt child = children[i]; 3209 PetscInt dof, off, j; 3210 3211 if (numSecFields) { 3212 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr); 3213 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr); 3214 } 3215 else { 3216 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr); 3217 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr); 3218 } 3219 3220 for (j = 0; j < dof / fComp; j++) { 3221 matCols[colOff++] = off + j * fComp; 3222 } 3223 } 3224 } 3225 if (classId == PETSCFE_CLASSID) { 3226 PetscFE fe = (PetscFE) disc; 3227 PetscInt fSize; 3228 3229 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 3230 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr); 3231 for (i = 0; i < numSelfShapes; i++) { /* for every shape function */ 3232 PetscQuadrature q; 3233 PetscInt dim, numPoints, j, k; 3234 const PetscReal *points; 3235 const PetscReal *weights; 3236 PetscInt *closure = NULL; 3237 PetscInt numClosure; 3238 PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfShapes - 1 - i) : i); 3239 PetscReal *Bparent; 3240 3241 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr); 3242 ierr = PetscQuadratureGetData(q,&dim,&numPoints,&points,&weights);CHKERRQ(ierr); 3243 ierr = PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3244 for (k = 0; k < numChildShapes; k++) { 3245 pointMat[numChildShapes * i + k] = 0.; 3246 } 3247 for (j = 0; j < numPoints; j++) { 3248 PetscInt childCell = -1; 3249 PetscReal parentValAtPoint; 3250 const PetscReal *pointReal = &points[dim * j]; 3251 const PetscScalar *point; 3252 PetscReal *Bchild; 3253 PetscInt childCellShapeOff, pointMatOff; 3254 #if defined(PETSC_USE_COMPLEX) 3255 PetscInt d; 3256 3257 for (d = 0; d < dim; d++) { 3258 pointScalar[d] = points[dim * j + d]; 3259 } 3260 point = pointScalar; 3261 #else 3262 point = pointReal; 3263 #endif 3264 3265 parentValAtPoint = Bparent[(fSize * j + parentCellShapeDof) * fComp]; 3266 3267 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3268 PetscInt child = children[k]; 3269 PetscInt *star = NULL; 3270 PetscInt numStar, s; 3271 3272 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3273 for (s = numStar - 1; s >= 0; s--) { 3274 PetscInt c = star[2 * s]; 3275 3276 if (c < cStart || c >= cEnd) continue; 3277 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr); 3278 if (childCell >= 0) break; 3279 } 3280 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3281 if (childCell >= 0) break; 3282 } 3283 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");CHKERRQ(ierr); 3284 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 3285 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr); 3286 CoordinatesRefToReal(dim, dim, v0parent, Jparent, pointReal, vtmp); 3287 CoordinatesRealToRef(dim, dim, v0, invJ, vtmp, pointRef); 3288 3289 ierr = PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3290 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3291 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3292 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3293 PetscInt l; 3294 3295 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr); 3296 childDof = depthNumDof[childDepth]; 3297 for (l = 0, childCellShapeOff = 0; l < numClosure; l++) { 3298 PetscInt point = closure[2 * l]; 3299 PetscInt pointDepth; 3300 3301 childO = closure[2 * l + 1]; 3302 if (point == child) break; 3303 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3304 childCellShapeOff += depthNumDof[pointDepth]; 3305 } 3306 if (l == numClosure) { 3307 pointMatOff += childDof; 3308 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3309 } 3310 for (l = 0; l < childDof; l++) { 3311 PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l); 3312 PetscReal childValAtPoint = Bchild[childCellDof * fComp]; 3313 3314 pointMat[i * numChildShapes + pointMatOff + l] += weights[j] * parentValAtPoint * childValAtPoint; 3315 } 3316 pointMatOff += childDof; 3317 } 3318 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3319 ierr = PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3320 } 3321 ierr = PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); 3322 } 3323 } 3324 else { /* just the volume-weighted averages of the children */ 3325 PetscInt childShapeOff; 3326 PetscReal parentVol; 3327 3328 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr); 3329 for (i = 0, childShapeOff = 0; i < numChildren; i++) { 3330 PetscInt child = children[i]; 3331 PetscReal childVol; 3332 3333 if (child < cStart || child >= cEnd) continue; 3334 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr); 3335 pointMat[childShapeOff] = childVol / parentVol; 3336 childShapeOff++; 3337 } 3338 } 3339 /* Insert pointMat into mat */ 3340 for (i = 0; i < fComp; i++) { 3341 PetscInt j; 3342 ierr = MatSetValues(mat,numSelfShapes,matRows,numChildShapes,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr); 3343 3344 for (j = 0; j < numSelfShapes; j++) { 3345 matRows[j]++; 3346 } 3347 for (j = 0; j < numChildShapes; j++) { 3348 matCols[j]++; 3349 } 3350 } 3351 ierr = DMRestoreWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3352 ierr = DMRestoreWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3353 } 3354 } 3355 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr); 3356 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr); 3357 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3358 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3359 *inj = mat; 3360 PetscFunctionReturn(0); 3361 } 3362 3363 #undef __FUNCT__ 3364 #define __FUNCT__ "DMPlexReferenceTreeGetChildrenMatrices_Injection" 3365 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3366 { 3367 PetscDS ds; 3368 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3369 PetscScalar ***refPointFieldMats; 3370 PetscSection refConSec, refSection; 3371 PetscErrorCode ierr; 3372 3373 PetscFunctionBegin; 3374 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3375 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3376 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3377 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3378 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3379 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 3380 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 3381 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 3382 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr); 3383 for (p = pRefStart; p < pRefEnd; p++) { 3384 PetscInt parent, pDof, parentDof; 3385 3386 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3387 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3388 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3389 if (!pDof || !parentDof || parent == p) continue; 3390 3391 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 3392 for (f = 0; f < numFields; f++) { 3393 PetscInt cDof, cOff, numCols, r, fComp; 3394 PetscObject disc; 3395 PetscClassId id; 3396 PetscFE fe = NULL; 3397 PetscFV fv = NULL; 3398 3399 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3400 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 3401 if (id == PETSCFE_CLASSID) { 3402 fe = (PetscFE) disc; 3403 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 3404 } 3405 else if (id == PETSCFV_CLASSID) { 3406 fv = (PetscFV) disc; 3407 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 3408 } 3409 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 3410 3411 if (numFields > 1) { 3412 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3413 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 3414 } 3415 else { 3416 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3417 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 3418 } 3419 3420 for (r = 0; r < cDof; r++) { 3421 rows[r] = cOff + r; 3422 } 3423 numCols = 0; 3424 { 3425 PetscInt aDof, aOff, j; 3426 3427 if (numFields > 1) { 3428 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr); 3429 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr); 3430 } 3431 else { 3432 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr); 3433 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr); 3434 } 3435 3436 for (j = 0; j < aDof; j++) { 3437 cols[numCols++] = aOff + j; 3438 } 3439 } 3440 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3441 /* transpose of constraint matrix */ 3442 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3443 } 3444 } 3445 *childrenMats = refPointFieldMats; 3446 ierr = PetscFree(rows);CHKERRQ(ierr); 3447 ierr = PetscFree(cols);CHKERRQ(ierr); 3448 PetscFunctionReturn(0); 3449 } 3450 3451 #undef __FUNCT__ 3452 #define __FUNCT__ "DMPlexReferenceTreeRestoreChildrenMatrices_Injection" 3453 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3454 { 3455 PetscDS ds; 3456 PetscScalar ***refPointFieldMats; 3457 PetscInt numFields, pRefStart, pRefEnd, p, f; 3458 PetscSection refConSec, refSection; 3459 PetscErrorCode ierr; 3460 3461 PetscFunctionBegin; 3462 refPointFieldMats = *childrenMats; 3463 *childrenMats = NULL; 3464 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3465 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3466 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3467 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3468 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3469 for (p = pRefStart; p < pRefEnd; p++) { 3470 PetscInt parent, pDof, parentDof; 3471 3472 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3473 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3474 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3475 if (!pDof || !parentDof || parent == p) continue; 3476 3477 for (f = 0; f < numFields; f++) { 3478 PetscInt cDof; 3479 3480 if (numFields > 1) { 3481 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3482 } 3483 else { 3484 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3485 } 3486 3487 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 3488 } 3489 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 3490 } 3491 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 3492 PetscFunctionReturn(0); 3493 } 3494 3495 #undef __FUNCT__ 3496 #define __FUNCT__ "DMPlexReferenceTreeGetInjector" 3497 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3498 { 3499 Mat cMatRef; 3500 PetscObject injRefObj; 3501 PetscErrorCode ierr; 3502 3503 PetscFunctionBegin; 3504 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr); 3505 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr); 3506 *injRef = (Mat) injRefObj; 3507 if (!*injRef) { 3508 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr); 3509 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr); 3510 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3511 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr); 3512 } 3513 PetscFunctionReturn(0); 3514 } 3515 3516 #undef __FUNCT__ 3517 #define __FUNCT__ "DMPlexTransferInjectorTree" 3518 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) 3519 { 3520 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3521 PetscSection globalCoarse, globalFine; 3522 PetscSection localCoarse, localFine, leafIndicesSec; 3523 PetscSection multiRootSec, rootIndicesSec; 3524 PetscInt *leafInds, *rootInds = NULL; 3525 const PetscInt *rootDegrees; 3526 PetscScalar *leafVals = NULL, *rootVals = NULL; 3527 PetscSF coarseToFineEmbedded; 3528 PetscErrorCode ierr; 3529 3530 PetscFunctionBegin; 3531 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3532 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3533 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3534 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3535 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 3536 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr); 3537 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr); 3538 { /* winnow fine points that don't have global dofs out of the sf */ 3539 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3540 const PetscInt *leaves; 3541 3542 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3543 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3544 p = leaves ? leaves[l] : l; 3545 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3546 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3547 if ((dof - cdof) > 0) { 3548 numPointsWithDofs++; 3549 3550 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr); 3551 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr); 3552 } 3553 } 3554 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3555 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr); 3556 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr); 3557 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr); 3558 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);} 3559 for (l = 0, offset = 0; l < nleaves; l++) { 3560 p = leaves ? leaves[l] : l; 3561 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3562 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3563 if ((dof - cdof) > 0) { 3564 PetscInt off, gOff; 3565 PetscInt *pInd; 3566 PetscScalar *pVal = NULL; 3567 3568 pointsWithDofs[offset++] = l; 3569 3570 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3571 3572 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3573 if (gatheredValues) { 3574 PetscInt i; 3575 3576 pVal = &leafVals[off + 1]; 3577 for (i = 0; i < dof; i++) pVal[i] = 0.; 3578 } 3579 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 3580 3581 offsets[0] = 0; 3582 if (numFields) { 3583 PetscInt f; 3584 3585 for (f = 0; f < numFields; f++) { 3586 PetscInt fDof; 3587 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr); 3588 offsets[f + 1] = fDof + offsets[f]; 3589 } 3590 DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 3591 } 3592 else { 3593 DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 3594 } 3595 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);} 3596 } 3597 } 3598 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 3599 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3600 } 3601 3602 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3603 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3604 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3605 3606 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3607 MPI_Datatype threeInt; 3608 PetscMPIInt rank; 3609 PetscInt (*parentNodeAndIdCoarse)[3]; 3610 PetscInt (*parentNodeAndIdFine)[3]; 3611 PetscInt p, nleaves, nleavesToParents; 3612 PetscSF pointSF, sfToParents; 3613 const PetscInt *ilocal; 3614 const PetscSFNode *iremote; 3615 PetscSFNode *iremoteToParents; 3616 PetscInt *ilocalToParents; 3617 3618 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr); 3619 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr); 3620 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr); 3621 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr); 3622 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr); 3623 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 3624 for (p = pStartC; p < pEndC; p++) { 3625 PetscInt parent, childId; 3626 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr); 3627 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3628 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3629 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3630 if (nleaves > 0) { 3631 PetscInt leaf = -1; 3632 3633 if (ilocal) { 3634 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr); 3635 } 3636 else { 3637 leaf = p - pStartC; 3638 } 3639 if (leaf >= 0) { 3640 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3641 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3642 } 3643 } 3644 } 3645 for (p = pStartF; p < pEndF; p++) { 3646 parentNodeAndIdFine[p - pStartF][0] = -1; 3647 parentNodeAndIdFine[p - pStartF][1] = -1; 3648 parentNodeAndIdFine[p - pStartF][2] = -1; 3649 } 3650 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3651 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3652 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3653 PetscInt dof; 3654 3655 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr); 3656 if (dof) { 3657 PetscInt off; 3658 3659 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3660 if (gatheredIndices) { 3661 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3662 } else if (gatheredValues) { 3663 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3664 } 3665 } 3666 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3667 nleavesToParents++; 3668 } 3669 } 3670 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr); 3671 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr); 3672 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3673 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3674 ilocalToParents[nleavesToParents] = p - pStartF; 3675 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3676 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3677 nleavesToParents++; 3678 } 3679 } 3680 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr); 3681 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr); 3682 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3683 3684 coarseToFineEmbedded = sfToParents; 3685 3686 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3687 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr); 3688 } 3689 3690 { /* winnow out coarse points that don't have dofs */ 3691 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3692 PetscSF sfDofsOnly; 3693 3694 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3695 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3696 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3697 if ((dof - cdof) > 0) { 3698 numPointsWithDofs++; 3699 } 3700 } 3701 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3702 for (p = pStartC, offset = 0; p < pEndC; p++) { 3703 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3704 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3705 if ((dof - cdof) > 0) { 3706 pointsWithDofs[offset++] = p - pStartC; 3707 } 3708 } 3709 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr); 3710 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3711 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3712 coarseToFineEmbedded = sfDofsOnly; 3713 } 3714 3715 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3716 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3717 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3718 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr); 3719 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr); 3720 for (p = pStartC; p < pEndC; p++) { 3721 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr); 3722 } 3723 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr); 3724 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr); 3725 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 3726 { /* distribute the leaf section */ 3727 PetscSF multi, multiInv, indicesSF; 3728 PetscInt *remoteOffsets, numRootIndices; 3729 3730 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr); 3731 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr); 3732 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr); 3733 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr); 3734 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 3735 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr); 3736 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 3737 if (gatheredIndices) { 3738 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr); 3739 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3740 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3741 } 3742 if (gatheredValues) { 3743 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr); 3744 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3745 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3746 } 3747 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 3748 } 3749 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3750 ierr = PetscFree(leafInds);CHKERRQ(ierr); 3751 ierr = PetscFree(leafVals);CHKERRQ(ierr); 3752 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3753 *rootMultiSec = multiRootSec; 3754 *multiLeafSec = rootIndicesSec; 3755 if (gatheredIndices) *gatheredIndices = rootInds; 3756 if (gatheredValues) *gatheredValues = rootVals; 3757 PetscFunctionReturn(0); 3758 } 3759 3760 #undef __FUNCT__ 3761 #define __FUNCT__ "DMPlexComputeInjectorTree" 3762 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3763 { 3764 DM refTree; 3765 PetscSection multiRootSec, rootIndicesSec; 3766 PetscSection globalCoarse, globalFine; 3767 PetscSection localCoarse, localFine; 3768 PetscSection cSecRef; 3769 PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd; 3770 Mat injRef; 3771 PetscInt numFields, maxDof; 3772 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3773 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3774 PetscLayout rowMap, colMap; 3775 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3776 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3777 PetscErrorCode ierr; 3778 3779 PetscFunctionBegin; 3780 3781 /* get the templates for the fine-to-coarse injection from the reference tree */ 3782 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 3783 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 3784 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3785 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 3786 3787 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3788 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3789 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3790 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 3791 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3792 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3793 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3794 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 3795 { 3796 PetscInt maxFields = PetscMax(1,numFields) + 1; 3797 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 3798 } 3799 3800 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr); 3801 3802 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr); 3803 3804 /* count indices */ 3805 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 3806 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 3807 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 3808 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 3809 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 3810 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr); 3811 for (p = pStartC; p < pEndC; p++) { 3812 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3813 3814 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3815 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3816 if ((dof - cdof) <= 0) continue; 3817 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3818 3819 rowOffsets[0] = 0; 3820 offsetsCopy[0] = 0; 3821 if (numFields) { 3822 PetscInt f; 3823 3824 for (f = 0; f < numFields; f++) { 3825 PetscInt fDof; 3826 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3827 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3828 } 3829 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3830 } 3831 else { 3832 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3833 rowOffsets[1] = offsetsCopy[0]; 3834 } 3835 3836 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3837 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3838 leafEnd = leafStart + numLeaves; 3839 for (l = leafStart; l < leafEnd; l++) { 3840 PetscInt numIndices, childId, offset; 3841 const PetscInt *childIndices; 3842 3843 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3844 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3845 childId = rootIndices[offset++]; 3846 childIndices = &rootIndices[offset]; 3847 numIndices--; 3848 3849 if (childId == -1) { /* equivalent points: scatter */ 3850 PetscInt i; 3851 3852 for (i = 0; i < numIndices; i++) { 3853 PetscInt colIndex = childIndices[i]; 3854 PetscInt rowIndex = parentIndices[i]; 3855 if (rowIndex < 0) continue; 3856 if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3857 if (colIndex >= colStart && colIndex < colEnd) { 3858 nnzD[rowIndex - rowStart] = 1; 3859 } 3860 else { 3861 nnzO[rowIndex - rowStart] = 1; 3862 } 3863 } 3864 } 3865 else { 3866 PetscInt parentId, f, lim; 3867 3868 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3869 3870 lim = PetscMax(1,numFields); 3871 offsets[0] = 0; 3872 if (numFields) { 3873 PetscInt f; 3874 3875 for (f = 0; f < numFields; f++) { 3876 PetscInt fDof; 3877 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3878 3879 offsets[f + 1] = fDof + offsets[f]; 3880 } 3881 } 3882 else { 3883 PetscInt cDof; 3884 3885 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3886 offsets[1] = cDof; 3887 } 3888 for (f = 0; f < lim; f++) { 3889 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3890 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3891 PetscInt i, numD = 0, numO = 0; 3892 3893 for (i = childStart; i < childEnd; i++) { 3894 PetscInt colIndex = childIndices[i]; 3895 3896 if (colIndex < 0) continue; 3897 if (colIndex >= colStart && colIndex < colEnd) { 3898 numD++; 3899 } 3900 else { 3901 numO++; 3902 } 3903 } 3904 for (i = parentStart; i < parentEnd; i++) { 3905 PetscInt rowIndex = parentIndices[i]; 3906 3907 if (rowIndex < 0) continue; 3908 nnzD[rowIndex - rowStart] += numD; 3909 nnzO[rowIndex - rowStart] += numO; 3910 } 3911 } 3912 } 3913 } 3914 } 3915 /* preallocate */ 3916 ierr = MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);CHKERRQ(ierr); 3917 ierr = PetscFree2(nnzD,nnzO);CHKERRQ(ierr); 3918 /* insert values */ 3919 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3920 for (p = pStartC; p < pEndC; p++) { 3921 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3922 3923 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3924 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3925 if ((dof - cdof) <= 0) continue; 3926 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3927 3928 rowOffsets[0] = 0; 3929 offsetsCopy[0] = 0; 3930 if (numFields) { 3931 PetscInt f; 3932 3933 for (f = 0; f < numFields; f++) { 3934 PetscInt fDof; 3935 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3936 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3937 } 3938 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3939 } 3940 else { 3941 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3942 rowOffsets[1] = offsetsCopy[0]; 3943 } 3944 3945 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3946 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3947 leafEnd = leafStart + numLeaves; 3948 for (l = leafStart; l < leafEnd; l++) { 3949 PetscInt numIndices, childId, offset; 3950 const PetscInt *childIndices; 3951 3952 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3953 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3954 childId = rootIndices[offset++]; 3955 childIndices = &rootIndices[offset]; 3956 numIndices--; 3957 3958 if (childId == -1) { /* equivalent points: scatter */ 3959 PetscInt i; 3960 3961 for (i = 0; i < numIndices; i++) { 3962 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr); 3963 } 3964 } 3965 else { 3966 PetscInt parentId, f, lim; 3967 3968 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3969 3970 lim = PetscMax(1,numFields); 3971 offsets[0] = 0; 3972 if (numFields) { 3973 PetscInt f; 3974 3975 for (f = 0; f < numFields; f++) { 3976 PetscInt fDof; 3977 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3978 3979 offsets[f + 1] = fDof + offsets[f]; 3980 } 3981 } 3982 else { 3983 PetscInt cDof; 3984 3985 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3986 offsets[1] = cDof; 3987 } 3988 for (f = 0; f < lim; f++) { 3989 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3990 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3991 const PetscInt *colIndices = &childIndices[offsets[f]]; 3992 3993 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr); 3994 } 3995 } 3996 } 3997 } 3998 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 3999 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4000 ierr = PetscFree(parentIndices);CHKERRQ(ierr); 4001 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4002 ierr = PetscFree(rootIndices);CHKERRQ(ierr); 4003 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4004 4005 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4006 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4007 PetscFunctionReturn(0); 4008 } 4009 4010 #undef __FUNCT__ 4011 #define __FUNCT__ "DMPlexTransferVecTree_Interpolate" 4012 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 4013 { 4014 PetscDS ds; 4015 PetscSF coarseToFineEmbedded; 4016 PetscSection globalCoarse, globalFine; 4017 PetscSection localCoarse, localFine; 4018 PetscSection aSec, cSec; 4019 PetscSection rootValuesSec; 4020 PetscSection leafValuesSec; 4021 PetscScalar *rootValues, *leafValues; 4022 IS aIS; 4023 const PetscInt *anchors; 4024 Mat cMat; 4025 PetscInt numFields; 4026 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior; 4027 PetscInt aStart, aEnd, cStart, cEnd; 4028 PetscInt *maxChildIds; 4029 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 4030 PetscFV fv = NULL; 4031 PetscInt dim, numFVcomps = -1, fvField = -1; 4032 DM cellDM = NULL, gradDM = NULL; 4033 const PetscScalar *cellGeomArray = NULL; 4034 const PetscScalar *gradArray = NULL; 4035 PetscErrorCode ierr; 4036 4037 PetscFunctionBegin; 4038 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4039 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4040 ierr = DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4041 ierr = DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4042 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4043 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4044 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4045 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr); 4046 { /* winnow fine points that don't have global dofs out of the sf */ 4047 PetscInt nleaves, l; 4048 const PetscInt *leaves; 4049 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 4050 4051 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 4052 4053 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 4054 PetscInt p = leaves ? leaves[l] : l; 4055 4056 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4057 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4058 if ((dof - cdof) > 0) { 4059 numPointsWithDofs++; 4060 } 4061 } 4062 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 4063 for (l = 0, offset = 0; l < nleaves; l++) { 4064 PetscInt p = leaves ? leaves[l] : l; 4065 4066 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4067 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4068 if ((dof - cdof) > 0) { 4069 pointsWithDofs[offset++] = l; 4070 } 4071 } 4072 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 4073 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 4074 } 4075 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4076 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 4077 for (p = pStartC; p < pEndC; p++) { 4078 maxChildIds[p - pStartC] = -2; 4079 } 4080 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4081 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4082 4083 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4084 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4085 4086 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 4087 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 4088 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 4089 4090 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 4091 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 4092 4093 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4094 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr); 4095 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr); 4096 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 4097 { 4098 PetscInt maxFields = PetscMax(1,numFields) + 1; 4099 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 4100 } 4101 if (grad) { 4102 PetscInt i; 4103 4104 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr); 4105 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr); 4106 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr); 4107 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr); 4108 for (i = 0; i < PetscMax(1,numFields); i++) { 4109 PetscObject obj; 4110 PetscClassId id; 4111 4112 ierr = DMGetField(coarse, i, &obj);CHKERRQ(ierr); 4113 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr); 4114 if (id == PETSCFV_CLASSID) { 4115 fv = (PetscFV) obj; 4116 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr); 4117 fvField = i; 4118 break; 4119 } 4120 } 4121 } 4122 4123 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4124 PetscInt dof; 4125 PetscInt maxChildId = maxChildIds[p - pStartC]; 4126 PetscInt numValues = 0; 4127 4128 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4129 if (dof < 0) { 4130 dof = -(dof + 1); 4131 } 4132 offsets[0] = 0; 4133 newOffsets[0] = 0; 4134 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4135 PetscInt *closure = NULL, closureSize, cl; 4136 4137 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4138 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4139 PetscInt c = closure[2 * cl], clDof; 4140 4141 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 4142 numValues += clDof; 4143 } 4144 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4145 } 4146 else if (maxChildId == -1) { 4147 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr); 4148 } 4149 /* we will pack the column indices with the field offsets */ 4150 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4151 /* also send the centroid, and the gradient */ 4152 numValues += dim * (1 + numFVcomps); 4153 } 4154 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr); 4155 } 4156 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr); 4157 { 4158 PetscInt numRootValues; 4159 const PetscScalar *coarseArray; 4160 4161 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr); 4162 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr); 4163 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4164 for (p = pStartC; p < pEndC; p++) { 4165 PetscInt numValues; 4166 PetscInt pValOff; 4167 PetscScalar *pVal; 4168 PetscInt maxChildId = maxChildIds[p - pStartC]; 4169 4170 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr); 4171 if (!numValues) { 4172 continue; 4173 } 4174 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr); 4175 pVal = &(rootValues[pValOff]); 4176 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4177 PetscInt closureSize = numValues; 4178 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr); 4179 if (grad && p >= cellStart && p < cellEnd) { 4180 PetscFVCellGeom *cg; 4181 PetscScalar *gradVals = NULL; 4182 PetscInt i; 4183 4184 pVal += (numValues - dim * (1 + numFVcomps)); 4185 4186 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr); 4187 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4188 pVal += dim; 4189 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr); 4190 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4191 } 4192 } 4193 else if (maxChildId == -1) { 4194 PetscInt lDof, lOff, i; 4195 4196 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr); 4197 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr); 4198 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4199 } 4200 } 4201 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4202 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 4203 } 4204 { 4205 PetscSF valuesSF; 4206 PetscInt *remoteOffsetsValues, numLeafValues; 4207 4208 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr); 4209 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr); 4210 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr); 4211 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 4212 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr); 4213 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr); 4214 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr); 4215 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4216 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4217 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr); 4218 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4219 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr); 4220 } 4221 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4222 { 4223 PetscInt maxDof; 4224 PetscInt *rowIndices; 4225 DM refTree; 4226 PetscInt **refPointFieldN; 4227 PetscScalar ***refPointFieldMats; 4228 PetscSection refConSec, refAnSec; 4229 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4230 PetscScalar *pointWork; 4231 4232 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 4233 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4234 ierr = DMGetWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4235 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 4236 ierr = DMGetDS(fine,&ds);CHKERRQ(ierr); 4237 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4238 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4239 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 4240 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 4241 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4242 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 4243 ierr = DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4244 ierr = DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4245 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4246 for (p = leafStart; p < leafEnd; p++) { 4247 PetscInt gDof, gcDof, gOff, lDof; 4248 PetscInt numValues, pValOff; 4249 PetscInt childId; 4250 const PetscScalar *pVal; 4251 const PetscScalar *fvGradData = NULL; 4252 4253 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 4254 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr); 4255 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 4256 if ((gDof - gcDof) <= 0) { 4257 continue; 4258 } 4259 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 4260 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr); 4261 if (!numValues) continue; 4262 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr); 4263 pVal = &leafValues[pValOff]; 4264 offsets[0] = 0; 4265 offsetsCopy[0] = 0; 4266 newOffsets[0] = 0; 4267 newOffsetsCopy[0] = 0; 4268 childId = cids[p - pStartF]; 4269 if (numFields) { 4270 PetscInt f; 4271 for (f = 0; f < numFields; f++) { 4272 PetscInt rowDof; 4273 4274 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 4275 offsets[f + 1] = offsets[f] + rowDof; 4276 offsetsCopy[f + 1] = offsets[f + 1]; 4277 /* TODO: closure indices */ 4278 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4279 } 4280 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 4281 } 4282 else { 4283 offsets[0] = 0; 4284 offsets[1] = lDof; 4285 newOffsets[0] = 0; 4286 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4287 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 4288 } 4289 if (childId == -1) { /* no child interpolation: one nnz per */ 4290 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr); 4291 } else { 4292 PetscInt f; 4293 4294 if (grad && p >= cellStart && p < cellEnd) { 4295 numValues -= (dim * (1 + numFVcomps)); 4296 fvGradData = &pVal[numValues]; 4297 } 4298 for (f = 0; f < PetscMax(1,numFields); f++) { 4299 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4300 PetscInt numRows = offsets[f+1] - offsets[f]; 4301 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4302 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4303 PetscScalar *rVal = &pointWork[offsets[f]]; 4304 PetscInt i, j; 4305 4306 #if 0 4307 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); 4308 #endif 4309 for (i = 0; i < numRows; i++) { 4310 PetscScalar val = 0.; 4311 for (j = 0; j < numCols; j++) { 4312 val += childMat[i * numCols + j] * cVal[j]; 4313 } 4314 rVal[i] = val; 4315 } 4316 if (f == fvField && p >= cellStart && p < cellEnd) { 4317 PetscReal centroid[3]; 4318 PetscScalar diff[3]; 4319 const PetscScalar *parentCentroid = &fvGradData[0]; 4320 const PetscScalar *gradient = &fvGradData[dim]; 4321 4322 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr); 4323 for (i = 0; i < dim; i++) { 4324 diff[i] = centroid[i] - parentCentroid[i]; 4325 } 4326 for (i = 0; i < numFVcomps; i++) { 4327 PetscScalar val = 0.; 4328 4329 for (j = 0; j < dim; j++) { 4330 val += gradient[dim * i + j] * diff[j]; 4331 } 4332 rVal[i] += val; 4333 } 4334 } 4335 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr); 4336 } 4337 } 4338 } 4339 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4340 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4341 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4342 } 4343 ierr = PetscFree(leafValues);CHKERRQ(ierr); 4344 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr); 4345 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 4346 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 4347 PetscFunctionReturn(0); 4348 } 4349 4350 #undef __FUNCT__ 4351 #define __FUNCT__ "DMPlexTransferVecTree_Inject" 4352 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4353 { 4354 PetscDS ds; 4355 DM refTree; 4356 PetscSection multiRootSec, rootIndicesSec; 4357 PetscSection globalCoarse, globalFine; 4358 PetscSection localCoarse, localFine; 4359 PetscSection cSecRef; 4360 PetscInt *parentIndices, pRefStart, pRefEnd; 4361 PetscScalar *rootValues, *parentValues; 4362 Mat injRef; 4363 PetscInt numFields, maxDof; 4364 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4365 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4366 PetscLayout rowMap, colMap; 4367 PetscInt rowStart, rowEnd, colStart, colEnd; 4368 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4369 PetscErrorCode ierr; 4370 4371 PetscFunctionBegin; 4372 4373 /* get the templates for the fine-to-coarse injection from the reference tree */ 4374 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4375 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4376 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 4377 ierr = DMGetDS(coarse,&ds);CHKERRQ(ierr); 4378 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4379 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 4380 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4381 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 4382 4383 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4384 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4385 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4386 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 4387 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4388 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4389 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4390 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 4391 { 4392 PetscInt maxFields = PetscMax(1,numFields) + 1; 4393 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 4394 } 4395 4396 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr); 4397 4398 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr); 4399 4400 /* count indices */ 4401 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr); 4402 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr); 4403 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 4404 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 4405 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 4406 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 4407 /* insert values */ 4408 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4409 for (p = pStartC; p < pEndC; p++) { 4410 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4411 PetscBool contribute = PETSC_FALSE; 4412 4413 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4414 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 4415 if ((dof - cdof) <= 0) continue; 4416 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr); 4417 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 4418 4419 rowOffsets[0] = 0; 4420 offsetsCopy[0] = 0; 4421 if (numFields) { 4422 PetscInt f; 4423 4424 for (f = 0; f < numFields; f++) { 4425 PetscInt fDof; 4426 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 4427 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4428 } 4429 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 4430 } 4431 else { 4432 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 4433 rowOffsets[1] = offsetsCopy[0]; 4434 } 4435 4436 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 4437 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 4438 leafEnd = leafStart + numLeaves; 4439 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4440 for (l = leafStart; l < leafEnd; l++) { 4441 PetscInt numIndices, childId, offset; 4442 const PetscScalar *childValues; 4443 4444 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 4445 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 4446 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4447 childValues = &rootValues[offset]; 4448 numIndices--; 4449 4450 if (childId == -2) { /* skip */ 4451 continue; 4452 } else if (childId == -1) { /* equivalent points: scatter */ 4453 PetscInt m; 4454 4455 contribute = PETSC_TRUE; 4456 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4457 } else { /* contributions from children: sum with injectors from reference tree */ 4458 PetscInt parentId, f, lim; 4459 4460 contribute = PETSC_TRUE; 4461 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 4462 4463 lim = PetscMax(1,numFields); 4464 offsets[0] = 0; 4465 if (numFields) { 4466 PetscInt f; 4467 4468 for (f = 0; f < numFields; f++) { 4469 PetscInt fDof; 4470 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4471 4472 offsets[f + 1] = fDof + offsets[f]; 4473 } 4474 } 4475 else { 4476 PetscInt cDof; 4477 4478 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4479 offsets[1] = cDof; 4480 } 4481 for (f = 0; f < lim; f++) { 4482 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4483 PetscInt n = offsets[f+1]-offsets[f]; 4484 PetscInt i, j; 4485 const PetscScalar *colValues = &childValues[offsets[f]]; 4486 4487 for (i = rowOffsets[f]; i < rowOffsets[f + 1]; i++) { 4488 PetscScalar val = 0.; 4489 for (j = 0; j < n; j++) { 4490 val += childMat[n * i + j] * colValues[j]; 4491 } 4492 parentValues[i] += val; 4493 } 4494 } 4495 } 4496 } 4497 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);} 4498 } 4499 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4500 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4501 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr); 4502 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4503 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4504 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4505 PetscFunctionReturn(0); 4506 } 4507 4508 #undef __FUNCT__ 4509 #define __FUNCT__ "DMPlexTransferVecTree" 4510 /*@ 4511 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4512 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4513 coarsening and refinement at the same time. 4514 4515 collective 4516 4517 Input Parameters: 4518 + dmIn - The DMPlex mesh for the input vector 4519 . vecIn - The input vector 4520 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4521 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4522 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4523 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4524 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4525 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4526 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4527 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4528 point j in dmOut is not a leaf of sfRefine. 4529 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4530 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4531 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4532 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4533 - time - Used if boundary values are time dependent. 4534 4535 Output Parameters: 4536 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered 4537 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4538 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4539 coarse points to fine points. 4540 4541 Level: developer 4542 4543 .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients() 4544 @*/ 4545 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4546 { 4547 PetscErrorCode ierr; 4548 4549 PetscFunctionBegin; 4550 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4551 if (sfRefine) { 4552 Vec vecInLocal; 4553 DM dmGrad = NULL; 4554 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4555 4556 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4557 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4558 { 4559 PetscInt numFields, i; 4560 4561 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4562 for (i = 0; i < numFields; i++) { 4563 PetscObject obj; 4564 PetscClassId classid; 4565 4566 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4567 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4568 if (classid == PETSCFV_CLASSID) { 4569 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4570 break; 4571 } 4572 } 4573 } 4574 if (useBCs) { 4575 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4576 } 4577 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4578 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4579 if (dmGrad) { 4580 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4581 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4582 } 4583 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr); 4584 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4585 if (dmGrad) { 4586 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4587 } 4588 } 4589 if (sfCoarsen) { 4590 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr); 4591 } 4592 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4593 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4594 PetscFunctionReturn(0); 4595 } 4596