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 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__ "DMPlexCreateDefaultReferenceTree" 69 /*@ 70 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 71 72 Collective on comm 73 74 Input Parameters: 75 + comm - the MPI communicator 76 . dim - the spatial dimension 77 - simplex - Flag for simplex, otherwise use a tensor-product cell 78 79 Output Parameters: 80 . ref - the reference tree DMPlex object 81 82 Level: intermediate 83 84 .keywords: reference cell 85 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 86 @*/ 87 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 88 { 89 DM K, Kref; 90 PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 91 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 92 DMLabel identity, identityRef; 93 PetscSection unionSection, unionConeSection, parentSection; 94 PetscScalar *unionCoords; 95 IS perm; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 /* create a reference element */ 100 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 101 ierr = DMPlexCreateLabel(K, "identity");CHKERRQ(ierr); 102 ierr = DMPlexGetLabel(K, "identity", &identity);CHKERRQ(ierr); 103 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 104 for (p = pStart; p < pEnd; p++) { 105 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 106 } 107 /* refine it */ 108 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 109 110 /* the reference tree is the union of these two, without duplicating 111 * points that appear in both */ 112 ierr = DMPlexGetLabel(Kref, "identity", &identityRef);CHKERRQ(ierr); 113 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 114 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 115 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 116 /* count points that will go in the union */ 117 for (p = pStart; p < pEnd; p++) { 118 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 119 } 120 for (p = pRefStart; p < pRefEnd; p++) { 121 PetscInt q, qSize; 122 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 123 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 124 if (qSize > 1) { 125 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 126 } 127 } 128 ierr = PetscMalloc1((pEnd - pStart) + (pRefEnd - pRefStart),&permvals);CHKERRQ(ierr); 129 offset = 0; 130 /* stratify points in the union by topological dimension */ 131 for (d = 0; d <= dim; d++) { 132 PetscInt cStart, cEnd, c; 133 134 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 135 for (c = cStart; c < cEnd; c++) { 136 permvals[offset++] = c; 137 } 138 139 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 140 for (c = cStart; c < cEnd; c++) { 141 permvals[offset++] = c + (pEnd - pStart); 142 } 143 } 144 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 145 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 146 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 147 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 148 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 149 /* count dimension points */ 150 for (d = 0; d <= dim; d++) { 151 PetscInt cStart, cOff, cOff2; 152 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 153 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 154 if (d < dim) { 155 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 156 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 157 } 158 else { 159 cOff2 = numUnionPoints; 160 } 161 numDimPoints[dim - d] = cOff2 - cOff; 162 } 163 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 164 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 165 /* count the cones in the union */ 166 for (p = pStart; p < pEnd; p++) { 167 PetscInt dof, uOff; 168 169 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 170 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 171 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 172 coneSizes[uOff] = dof; 173 } 174 for (p = pRefStart; p < pRefEnd; p++) { 175 PetscInt dof, uDof, uOff; 176 177 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 178 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 179 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 180 if (uDof) { 181 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 182 coneSizes[uOff] = dof; 183 } 184 } 185 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 186 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 187 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 188 /* write the cones in the union */ 189 for (p = pStart; p < pEnd; p++) { 190 PetscInt dof, uOff, c, cOff; 191 const PetscInt *cone, *orientation; 192 193 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 194 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 195 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 196 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 197 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 198 for (c = 0; c < dof; c++) { 199 PetscInt e, eOff; 200 e = cone[c]; 201 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 202 unionCones[cOff + c] = eOff; 203 unionOrientations[cOff + c] = orientation[c]; 204 } 205 } 206 for (p = pRefStart; p < pRefEnd; p++) { 207 PetscInt dof, uDof, uOff, c, cOff; 208 const PetscInt *cone, *orientation; 209 210 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 211 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 212 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 213 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 214 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 215 if (uDof) { 216 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 217 for (c = 0; c < dof; c++) { 218 PetscInt e, eOff, eDof; 219 220 e = cone[c]; 221 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 222 if (eDof) { 223 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 224 } 225 else { 226 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 227 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 228 } 229 unionCones[cOff + c] = eOff; 230 unionOrientations[cOff + c] = orientation[c]; 231 } 232 } 233 } 234 /* get the coordinates */ 235 { 236 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 237 PetscSection KcoordsSec, KrefCoordsSec; 238 Vec KcoordsVec, KrefCoordsVec; 239 PetscScalar *Kcoords; 240 241 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 242 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 243 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 244 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 245 246 numVerts = numDimPoints[0]; 247 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 248 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 249 250 offset = 0; 251 for (v = vStart; v < vEnd; v++) { 252 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 253 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 254 for (d = 0; d < dim; d++) { 255 unionCoords[offset * dim + d] = Kcoords[d]; 256 } 257 offset++; 258 } 259 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 260 for (v = vRefStart; v < vRefEnd; v++) { 261 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 262 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 263 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 264 if (vDof) { 265 for (d = 0; d < dim; d++) { 266 unionCoords[offset * dim + d] = Kcoords[d]; 267 } 268 offset++; 269 } 270 } 271 } 272 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 273 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 274 ierr = DMPlexSetDimension(*ref,dim);CHKERRQ(ierr); 275 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 276 /* set the tree */ 277 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 278 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 279 for (p = pRefStart; p < pRefEnd; p++) { 280 PetscInt uDof, uOff; 281 282 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 283 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 284 if (uDof) { 285 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 286 } 287 } 288 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 289 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 290 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 291 for (p = pRefStart; p < pRefEnd; p++) { 292 PetscInt uDof, uOff; 293 294 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 295 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 296 if (uDof) { 297 PetscInt pOff, parent, parentU; 298 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 299 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 300 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 301 parents[pOff] = parentU; 302 childIDs[pOff] = uOff; 303 } 304 } 305 ierr = DMPlexSetTree(*ref,parentSection,parents,childIDs);CHKERRQ(ierr); 306 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 307 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 308 309 /* clean up */ 310 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 311 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 312 ierr = ISDestroy(&perm);CHKERRQ(ierr); 313 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 314 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 315 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 316 ierr = DMDestroy(&K);CHKERRQ(ierr); 317 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 } 320 321 #undef __FUNCT__ 322 #define __FUNCT__ "DMPlexTreeSymmetrize" 323 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 324 { 325 DM_Plex *mesh = (DM_Plex *)dm->data; 326 PetscSection childSec, pSec; 327 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 328 PetscInt *offsets, *children, pStart, pEnd; 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 333 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 334 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 335 pSec = mesh->parentSection; 336 if (!pSec) PetscFunctionReturn(0); 337 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 338 for (p = 0; p < pSize; p++) { 339 PetscInt par = mesh->parents[p]; 340 341 parMax = PetscMax(parMax,par+1); 342 parMin = PetscMin(parMin,par); 343 } 344 if (parMin > parMax) { 345 parMin = -1; 346 parMax = -1; 347 } 348 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 349 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 350 for (p = 0; p < pSize; p++) { 351 PetscInt par = mesh->parents[p]; 352 353 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 354 } 355 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 356 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 357 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 358 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 359 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 360 for (p = pStart; p < pEnd; p++) { 361 PetscInt dof, off, i; 362 363 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 364 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 365 for (i = 0; i < dof; i++) { 366 PetscInt par = mesh->parents[off + i], cOff; 367 368 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 369 children[cOff + offsets[par-parMin]++] = p; 370 } 371 } 372 mesh->childSection = childSec; 373 mesh->children = children; 374 ierr = PetscFree(offsets);CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "AnchorsFlatten" 380 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 381 { 382 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 383 const PetscInt *vals; 384 PetscSection secNew; 385 PetscBool anyNew, globalAnyNew; 386 PetscBool compress; 387 PetscErrorCode ierr; 388 389 PetscFunctionBegin; 390 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 391 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 392 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 393 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 394 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 395 for (i = 0; i < size; i++) { 396 PetscInt dof; 397 398 p = vals[i]; 399 if (p < pStart || p >= pEnd) continue; 400 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 401 if (dof) break; 402 } 403 if (i == size) { 404 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 405 anyNew = PETSC_FALSE; 406 compress = PETSC_FALSE; 407 sizeNew = 0; 408 } 409 else { 410 anyNew = PETSC_TRUE; 411 for (p = pStart; p < pEnd; p++) { 412 PetscInt dof, off; 413 414 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 415 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 416 for (i = 0; i < dof; i++) { 417 PetscInt q = vals[off + i], qDof = 0; 418 419 if (q >= pStart && q < pEnd) { 420 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 421 } 422 if (qDof) { 423 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 424 } 425 else { 426 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 427 } 428 } 429 } 430 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 431 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 432 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 433 compress = PETSC_FALSE; 434 for (p = pStart; p < pEnd; p++) { 435 PetscInt dof, off, count, offNew, dofNew; 436 437 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 438 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 439 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 440 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 441 count = 0; 442 for (i = 0; i < dof; i++) { 443 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 444 445 if (q >= pStart && q < pEnd) { 446 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 447 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 448 } 449 if (qDof) { 450 PetscInt oldCount = count; 451 452 for (j = 0; j < qDof; j++) { 453 PetscInt k, r = vals[qOff + j]; 454 455 for (k = 0; k < oldCount; k++) { 456 if (valsNew[offNew + k] == r) { 457 break; 458 } 459 } 460 if (k == oldCount) { 461 valsNew[offNew + count++] = r; 462 } 463 } 464 } 465 else { 466 PetscInt k, oldCount = count; 467 468 for (k = 0; k < oldCount; k++) { 469 if (valsNew[offNew + k] == q) { 470 break; 471 } 472 } 473 if (k == oldCount) { 474 valsNew[offNew + count++] = q; 475 } 476 } 477 } 478 if (count < dofNew) { 479 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 480 compress = PETSC_TRUE; 481 } 482 } 483 } 484 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 485 ierr = MPI_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 486 if (!globalAnyNew) { 487 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 488 *sectionNew = NULL; 489 *isNew = NULL; 490 } 491 else { 492 PetscBool globalCompress; 493 494 ierr = MPI_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 495 if (compress) { 496 PetscSection secComp; 497 PetscInt *valsComp = NULL; 498 499 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 500 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 501 for (p = pStart; p < pEnd; p++) { 502 PetscInt dof; 503 504 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 505 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 506 } 507 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 508 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 509 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 510 for (p = pStart; p < pEnd; p++) { 511 PetscInt dof, off, offNew, j; 512 513 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 514 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 515 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 516 for (j = 0; j < dof; j++) { 517 valsComp[offNew + j] = valsNew[off + j]; 518 } 519 } 520 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 521 secNew = secComp; 522 ierr = PetscFree(valsNew);CHKERRQ(ierr); 523 valsNew = valsComp; 524 } 525 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 526 } 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "DMPlexComputeConstraints_Tree" 532 static PetscErrorCode DMPlexComputeConstraints_Tree(DM dm) 533 { 534 PetscInt p, pStart, pEnd, *anchors, size; 535 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 536 PetscSection aSec; 537 IS aIS; 538 PetscErrorCode ierr; 539 540 PetscFunctionBegin; 541 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 542 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 543 for (p = pStart; p < pEnd; p++) { 544 PetscInt parent; 545 546 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 547 if (parent != p) { 548 aMin = PetscMin(aMin,p); 549 aMax = PetscMax(aMax,p+1); 550 } 551 } 552 if (aMin > aMax) { 553 aMin = -1; 554 aMax = -1; 555 } 556 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&aSec);CHKERRQ(ierr); 557 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 558 for (p = aMin; p < aMax; p++) { 559 PetscInt parent, ancestor = p; 560 561 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 562 while (parent != ancestor) { 563 ancestor = parent; 564 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 565 } 566 if (ancestor != p) { 567 PetscInt closureSize, *closure = NULL; 568 569 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 570 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 571 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 572 } 573 } 574 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 575 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 576 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 577 for (p = aMin; p < aMax; p++) { 578 PetscInt parent, ancestor = p; 579 580 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 581 while (parent != ancestor) { 582 ancestor = parent; 583 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 584 } 585 if (ancestor != p) { 586 PetscInt j, closureSize, *closure = NULL, aOff; 587 588 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 589 590 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 591 for (j = 0; j < closureSize; j++) { 592 anchors[aOff + j] = closure[2*j]; 593 } 594 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 595 } 596 } 597 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 598 { 599 PetscSection aSecNew = aSec; 600 IS aISNew = aIS; 601 602 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 603 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 604 while (aSecNew) { 605 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 606 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 607 aSec = aSecNew; 608 aIS = aISNew; 609 aSecNew = NULL; 610 aISNew = NULL; 611 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 612 } 613 } 614 ierr = DMPlexSetConstraints(dm,aSec,aIS);CHKERRQ(ierr); 615 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 616 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "DMPlexSetTree" 622 /*@ 623 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 624 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 625 tree root. 626 627 Collective on dm 628 629 Input Parameters: 630 + dm - the DMPlex object 631 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 632 offset indexes the parent and childID list; the reference count of parentSection is incremented 633 . parents - a list of the point parents; copied, can be destroyed 634 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 635 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 636 637 Level: intermediate 638 639 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 640 @*/ 641 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 642 { 643 DM_Plex *mesh = (DM_Plex *)dm->data; 644 PetscInt size; 645 PetscErrorCode ierr; 646 647 PetscFunctionBegin; 648 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 649 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 650 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 651 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 652 mesh->parentSection = parentSection; 653 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 654 if (parents != mesh->parents) { 655 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 656 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 657 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 658 } 659 if (childIDs != mesh->childIDs) { 660 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 661 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 662 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 663 } 664 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 665 ierr = DMPlexComputeConstraints_Tree(dm);CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "DMPlexGetTree" 671 /*@ 672 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 673 Collective on dm 674 675 Input Parameters: 676 . dm - the DMPlex object 677 678 Output Parameters: 679 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 680 offset indexes the parent and childID list 681 . parents - a list of the point parents 682 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 683 the child corresponds to the point in the reference tree with index childID 684 . childSection - the inverse of the parent section 685 - children - a list of the point children 686 687 Level: intermediate 688 689 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 690 @*/ 691 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 692 { 693 DM_Plex *mesh = (DM_Plex *)dm->data; 694 695 PetscFunctionBegin; 696 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 697 if (parentSection) *parentSection = mesh->parentSection; 698 if (parents) *parents = mesh->parents; 699 if (childIDs) *childIDs = mesh->childIDs; 700 if (childSection) *childSection = mesh->childSection; 701 if (children) *children = mesh->children; 702 PetscFunctionReturn(0); 703 } 704 705 #undef __FUNCT__ 706 #define __FUNCT__ "DMPlexGetTreeParent" 707 /*@ 708 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 709 710 Input Parameters: 711 + dm - the DMPlex object 712 - point - the query point 713 714 Output Parameters: 715 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 716 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 717 does not have a parent 718 719 Level: intermediate 720 721 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 722 @*/ 723 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 724 { 725 DM_Plex *mesh = (DM_Plex *)dm->data; 726 PetscSection pSec; 727 PetscErrorCode ierr; 728 729 PetscFunctionBegin; 730 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 731 pSec = mesh->parentSection; 732 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 733 PetscInt dof; 734 735 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 736 if (dof) { 737 PetscInt off; 738 739 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 740 if (parent) *parent = mesh->parents[off]; 741 if (childID) *childID = mesh->childIDs[off]; 742 PetscFunctionReturn(0); 743 } 744 } 745 if (parent) { 746 *parent = point; 747 } 748 if (childID) { 749 *childID = 0; 750 } 751 PetscFunctionReturn(0); 752 } 753 754 #undef __FUNCT__ 755 #define __FUNCT__ "DMPlexGetTreeChildren" 756 /*@C 757 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 758 759 Input Parameters: 760 + dm - the DMPlex object 761 - point - the query point 762 763 Output Parameters: 764 + numChildren - if not NULL, set to the number of children 765 - children - if not NULL, set to a list children, or set to NULL if the point has no children 766 767 Level: intermediate 768 769 Fortran Notes: 770 Since it returns an array, this routine is only available in Fortran 90, and you must 771 include petsc.h90 in your code. 772 773 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 774 @*/ 775 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 776 { 777 DM_Plex *mesh = (DM_Plex *)dm->data; 778 PetscSection childSec; 779 PetscInt dof = 0; 780 PetscErrorCode ierr; 781 782 PetscFunctionBegin; 783 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 784 childSec = mesh->childSection; 785 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 786 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 787 } 788 if (numChildren) *numChildren = dof; 789 if (children) { 790 if (dof) { 791 PetscInt off; 792 793 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 794 *children = &mesh->children[off]; 795 } 796 else { 797 *children = NULL; 798 } 799 } 800 PetscFunctionReturn(0); 801 } 802 803 #undef __FUNCT__ 804 #define __FUNCT__ "DMPlexComputeConstraintMatrix_ReferenceTree" 805 PetscErrorCode DMPlexComputeConstraintMatrix_ReferenceTree(DM dm) 806 { 807 PetscDS ds; 808 PetscInt spdim; 809 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 810 const PetscInt *anchors; 811 PetscSection section, cSec, aSec; 812 Mat cMat; 813 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 814 IS aIS; 815 PetscErrorCode ierr; 816 817 PetscFunctionBegin; 818 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 819 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 820 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 821 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 822 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 823 ierr = DMPlexGetConstraintSection(dm,&cSec);CHKERRQ(ierr); 824 ierr = DMPlexGetConstraints(dm,&aSec,&aIS);CHKERRQ(ierr); 825 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 826 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 827 ierr = DMPlexGetConstraintMatrix(dm,&cMat);CHKERRQ(ierr); 828 ierr = DMPlexGetDimension(dm,&spdim);CHKERRQ(ierr); 829 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 830 831 for (f = 0; f < numFields; f++) { 832 PetscFE fe; 833 PetscDualSpace space; 834 PetscInt i, j, k, nPoints, offset; 835 PetscInt fSize, fComp; 836 PetscScalar *B = NULL; 837 PetscReal *weights, *pointsRef, *pointsReal; 838 Mat Amat, Bmat, Xmat; 839 840 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *)(&fe));CHKERRQ(ierr); 841 ierr = PetscFEGetDualSpace(fe,&space);CHKERRQ(ierr); 842 ierr = PetscDualSpaceGetDimension(space,&fSize);CHKERRQ(ierr); 843 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 844 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 845 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 846 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 847 ierr = MatSetUp(Amat);CHKERRQ(ierr); 848 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 849 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 850 nPoints = 0; 851 for (i = 0; i < fSize; i++) { 852 PetscInt qPoints; 853 PetscQuadrature quad; 854 855 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 856 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 857 nPoints += qPoints; 858 } 859 ierr = PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);CHKERRQ(ierr); 860 offset = 0; 861 for (i = 0; i < fSize; i++) { 862 PetscInt qPoints; 863 const PetscReal *p, *w; 864 PetscQuadrature quad; 865 866 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 867 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 868 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 869 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 870 offset += qPoints; 871 } 872 ierr = PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 873 offset = 0; 874 for (i = 0; i < fSize; i++) { 875 PetscInt qPoints; 876 PetscQuadrature quad; 877 878 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 879 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 880 for (j = 0; j < fSize; j++) { 881 PetscScalar val = 0.; 882 883 for (k = 0; k < qPoints; k++) { 884 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 885 } 886 ierr = MatSetValue(Amat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 887 } 888 offset += qPoints; 889 } 890 ierr = PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 891 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 892 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 893 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 894 for (c = cStart; c < cEnd; c++) { 895 PetscInt parent; 896 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 897 PetscInt *childOffsets, *parentOffsets; 898 899 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 900 if (parent == c) continue; 901 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 902 for (i = 0; i < closureSize; i++) { 903 PetscInt p = closure[2*i]; 904 PetscInt conDof; 905 906 if (p < conStart || p >= conEnd) continue; 907 if (numFields > 1) { 908 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 909 } 910 else { 911 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 912 } 913 if (conDof) break; 914 } 915 if (i == closureSize) { 916 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 917 continue; 918 } 919 920 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 921 ierr = DMPlexComputeCellGeometry(dm, parent, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 922 for (i = 0; i < nPoints; i++) { 923 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);CHKERRQ(ierr); 924 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);CHKERRQ(ierr); 925 } 926 ierr = PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 927 offset = 0; 928 for (i = 0; i < fSize; i++) { 929 PetscInt qPoints; 930 PetscQuadrature quad; 931 932 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 933 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 934 for (j = 0; j < fSize; j++) { 935 PetscScalar val = 0.; 936 937 for (k = 0; k < qPoints; k++) { 938 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 939 } 940 MatSetValue(Bmat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 941 } 942 offset += qPoints; 943 } 944 ierr = PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 945 ierr = MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 946 ierr = MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 947 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 948 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 949 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 950 childOffsets[0] = 0; 951 for (i = 0; i < closureSize; i++) { 952 PetscInt p = closure[2*i]; 953 PetscInt dof; 954 955 if (numFields > 1) { 956 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 957 } 958 else { 959 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 960 } 961 childOffsets[i+1]=childOffsets[i]+dof / fComp; 962 } 963 parentOffsets[0] = 0; 964 for (i = 0; i < closureSizeP; i++) { 965 PetscInt p = closureP[2*i]; 966 PetscInt dof; 967 968 if (numFields > 1) { 969 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 970 } 971 else { 972 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 973 } 974 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 975 } 976 for (i = 0; i < closureSize; i++) { 977 PetscInt conDof, conOff, aDof, aOff; 978 PetscInt p = closure[2*i]; 979 PetscInt o = closure[2*i+1]; 980 981 if (p < conStart || p >= conEnd) continue; 982 if (numFields > 1) { 983 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 984 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 985 } 986 else { 987 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 988 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 989 } 990 if (!conDof) continue; 991 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 992 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 993 for (k = 0; k < aDof; k++) { 994 PetscInt a = anchors[aOff + k]; 995 PetscInt aSecDof, aSecOff; 996 997 if (numFields > 1) { 998 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 999 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1000 } 1001 else { 1002 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1003 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1004 } 1005 if (!aSecDof) continue; 1006 1007 for (j = 0; j < closureSizeP; j++) { 1008 PetscInt q = closureP[2*j]; 1009 PetscInt oq = closureP[2*j+1]; 1010 1011 if (q == a) { 1012 PetscInt r, s, t; 1013 1014 for (r = childOffsets[i]; r < childOffsets[i+1]; r++) { 1015 for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) { 1016 PetscScalar val; 1017 PetscInt insertCol, insertRow; 1018 1019 ierr = MatGetValue(Xmat,r,s,&val);CHKERRQ(ierr); 1020 if (o >= 0) { 1021 insertRow = conOff + fComp * (r - childOffsets[i]); 1022 } 1023 else { 1024 insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r); 1025 } 1026 if (oq >= 0) { 1027 insertCol = aSecOff + fComp * (s - parentOffsets[j]); 1028 } 1029 else { 1030 insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s); 1031 } 1032 for (t = 0; t < fComp; t++) { 1033 ierr = MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);CHKERRQ(ierr); 1034 } 1035 } 1036 } 1037 } 1038 } 1039 } 1040 } 1041 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1042 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1043 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1044 } 1045 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1046 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1047 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1048 ierr = PetscFree3(weights,pointsRef,pointsReal);CHKERRQ(ierr); 1049 } 1050 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1051 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1052 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1053 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1054 1055 PetscFunctionReturn(0); 1056 } 1057