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 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfIn, PetscSF sfOut, PetscInt *cidsIn, PetscInt *cidsOut, PetscBool useBCs, PetscReal time) 4475 { 4476 PetscErrorCode ierr; 4477 4478 PetscFunctionBegin; 4479 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4480 if (sfIn) { 4481 Vec vecInLocal; 4482 DM dmGrad = NULL; 4483 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4484 4485 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4486 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4487 { 4488 PetscInt numFields, i; 4489 4490 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4491 for (i = 0; i < numFields; i++) { 4492 PetscObject obj; 4493 PetscClassId classid; 4494 4495 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4496 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4497 if (classid == PETSCFV_CLASSID) { 4498 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4499 break; 4500 } 4501 } 4502 } 4503 if (useBCs) { 4504 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4505 } 4506 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4507 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4508 if (dmGrad) { 4509 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4510 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4511 } 4512 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfIn,cidsIn,grad,cellGeom);CHKERRQ(ierr); 4513 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4514 if (dmGrad) { 4515 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4516 } 4517 } 4518 if (sfOut) { 4519 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfOut,cidsOut);CHKERRQ(ierr); 4520 } 4521 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4522 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4523 PetscFunctionReturn(0); 4524 } 4525