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