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