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 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool); 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "DMPlexCreateDefaultReferenceTree" 71 /*@ 72 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 73 74 Collective on comm 75 76 Input Parameters: 77 + comm - the MPI communicator 78 . dim - the spatial dimension 79 - simplex - Flag for simplex, otherwise use a tensor-product cell 80 81 Output Parameters: 82 . ref - the reference tree DMPlex object 83 84 Level: intermediate 85 86 .keywords: reference cell 87 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 88 @*/ 89 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 90 { 91 DM K, Kref; 92 PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 93 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 94 DMLabel identity, identityRef; 95 PetscSection unionSection, unionConeSection, parentSection; 96 PetscScalar *unionCoords; 97 IS perm; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 /* create a reference element */ 102 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 103 ierr = DMPlexCreateLabel(K, "identity");CHKERRQ(ierr); 104 ierr = DMPlexGetLabel(K, "identity", &identity);CHKERRQ(ierr); 105 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 106 for (p = pStart; p < pEnd; p++) { 107 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 108 } 109 /* refine it */ 110 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 111 112 /* the reference tree is the union of these two, without duplicating 113 * points that appear in both */ 114 ierr = DMPlexGetLabel(Kref, "identity", &identityRef);CHKERRQ(ierr); 115 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 116 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 117 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 118 /* count points that will go in the union */ 119 for (p = pStart; p < pEnd; p++) { 120 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 121 } 122 for (p = pRefStart; p < pRefEnd; p++) { 123 PetscInt q, qSize; 124 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 125 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 126 if (qSize > 1) { 127 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 128 } 129 } 130 ierr = PetscMalloc1((pEnd - pStart) + (pRefEnd - pRefStart),&permvals);CHKERRQ(ierr); 131 offset = 0; 132 /* stratify points in the union by topological dimension */ 133 for (d = 0; d <= dim; d++) { 134 PetscInt cStart, cEnd, c; 135 136 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 137 for (c = cStart; c < cEnd; c++) { 138 permvals[offset++] = c; 139 } 140 141 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 142 for (c = cStart; c < cEnd; c++) { 143 permvals[offset++] = c + (pEnd - pStart); 144 } 145 } 146 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 147 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 148 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 149 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 150 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 151 /* count dimension points */ 152 for (d = 0; d <= dim; d++) { 153 PetscInt cStart, cOff, cOff2; 154 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 155 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 156 if (d < dim) { 157 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 158 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 159 } 160 else { 161 cOff2 = numUnionPoints; 162 } 163 numDimPoints[dim - d] = cOff2 - cOff; 164 } 165 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 166 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 167 /* count the cones in the union */ 168 for (p = pStart; p < pEnd; p++) { 169 PetscInt dof, uOff; 170 171 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 172 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 173 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 174 coneSizes[uOff] = dof; 175 } 176 for (p = pRefStart; p < pRefEnd; p++) { 177 PetscInt dof, uDof, uOff; 178 179 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 180 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 181 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 182 if (uDof) { 183 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 184 coneSizes[uOff] = dof; 185 } 186 } 187 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 188 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 189 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 190 /* write the cones in the union */ 191 for (p = pStart; p < pEnd; p++) { 192 PetscInt dof, uOff, c, cOff; 193 const PetscInt *cone, *orientation; 194 195 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 196 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 197 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 198 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 199 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 200 for (c = 0; c < dof; c++) { 201 PetscInt e, eOff; 202 e = cone[c]; 203 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 204 unionCones[cOff + c] = eOff; 205 unionOrientations[cOff + c] = orientation[c]; 206 } 207 } 208 for (p = pRefStart; p < pRefEnd; p++) { 209 PetscInt dof, uDof, uOff, c, cOff; 210 const PetscInt *cone, *orientation; 211 212 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 213 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 214 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 215 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 216 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 217 if (uDof) { 218 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 219 for (c = 0; c < dof; c++) { 220 PetscInt e, eOff, eDof; 221 222 e = cone[c]; 223 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 224 if (eDof) { 225 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 226 } 227 else { 228 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 229 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 230 } 231 unionCones[cOff + c] = eOff; 232 unionOrientations[cOff + c] = orientation[c]; 233 } 234 } 235 } 236 /* get the coordinates */ 237 { 238 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 239 PetscSection KcoordsSec, KrefCoordsSec; 240 Vec KcoordsVec, KrefCoordsVec; 241 PetscScalar *Kcoords; 242 243 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 244 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 245 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 246 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 247 248 numVerts = numDimPoints[0]; 249 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 250 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 251 252 offset = 0; 253 for (v = vStart; v < vEnd; v++) { 254 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 255 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 256 for (d = 0; d < dim; d++) { 257 unionCoords[offset * dim + d] = Kcoords[d]; 258 } 259 offset++; 260 } 261 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 262 for (v = vRefStart; v < vRefEnd; v++) { 263 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 264 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 265 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 266 if (vDof) { 267 for (d = 0; d < dim; d++) { 268 unionCoords[offset * dim + d] = Kcoords[d]; 269 } 270 offset++; 271 } 272 } 273 } 274 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 275 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 276 ierr = DMPlexSetDimension(*ref,dim);CHKERRQ(ierr); 277 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 278 /* set the tree */ 279 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 280 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 281 for (p = pRefStart; p < pRefEnd; p++) { 282 PetscInt uDof, uOff; 283 284 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 285 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 286 if (uDof) { 287 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 288 } 289 } 290 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 291 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 292 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 293 for (p = pRefStart; p < pRefEnd; p++) { 294 PetscInt uDof, uOff; 295 296 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 297 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 298 if (uDof) { 299 PetscInt pOff, parent, parentU; 300 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 301 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 302 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 303 parents[pOff] = parentU; 304 childIDs[pOff] = uOff; 305 } 306 } 307 ierr = DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE);CHKERRQ(ierr); 308 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 309 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 310 311 /* clean up */ 312 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 313 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 314 ierr = ISDestroy(&perm);CHKERRQ(ierr); 315 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 316 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 317 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 318 ierr = DMDestroy(&K);CHKERRQ(ierr); 319 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "DMPlexTreeSymmetrize" 326 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 327 { 328 DM_Plex *mesh = (DM_Plex *)dm->data; 329 PetscSection childSec, pSec; 330 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 331 PetscInt *offsets, *children, pStart, pEnd; 332 PetscErrorCode ierr; 333 334 PetscFunctionBegin; 335 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 336 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 337 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 338 pSec = mesh->parentSection; 339 if (!pSec) PetscFunctionReturn(0); 340 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 341 for (p = 0; p < pSize; p++) { 342 PetscInt par = mesh->parents[p]; 343 344 parMax = PetscMax(parMax,par+1); 345 parMin = PetscMin(parMin,par); 346 } 347 if (parMin > parMax) { 348 parMin = -1; 349 parMax = -1; 350 } 351 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 352 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 353 for (p = 0; p < pSize; p++) { 354 PetscInt par = mesh->parents[p]; 355 356 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 357 } 358 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 359 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 360 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 361 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 362 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 363 for (p = pStart; p < pEnd; p++) { 364 PetscInt dof, off, i; 365 366 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 367 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 368 for (i = 0; i < dof; i++) { 369 PetscInt par = mesh->parents[off + i], cOff; 370 371 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 372 children[cOff + offsets[par-parMin]++] = p; 373 } 374 } 375 mesh->childSection = childSec; 376 mesh->children = children; 377 ierr = PetscFree(offsets);CHKERRQ(ierr); 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNCT__ 382 #define __FUNCT__ "AnchorsFlatten" 383 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 384 { 385 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 386 const PetscInt *vals; 387 PetscSection secNew; 388 PetscBool anyNew, globalAnyNew; 389 PetscBool compress; 390 PetscErrorCode ierr; 391 392 PetscFunctionBegin; 393 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 394 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 395 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 396 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 397 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 398 for (i = 0; i < size; i++) { 399 PetscInt dof; 400 401 p = vals[i]; 402 if (p < pStart || p >= pEnd) continue; 403 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 404 if (dof) break; 405 } 406 if (i == size) { 407 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 408 anyNew = PETSC_FALSE; 409 compress = PETSC_FALSE; 410 sizeNew = 0; 411 } 412 else { 413 anyNew = PETSC_TRUE; 414 for (p = pStart; p < pEnd; p++) { 415 PetscInt dof, off; 416 417 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 418 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 419 for (i = 0; i < dof; i++) { 420 PetscInt q = vals[off + i], qDof = 0; 421 422 if (q >= pStart && q < pEnd) { 423 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 424 } 425 if (qDof) { 426 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 427 } 428 else { 429 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 430 } 431 } 432 } 433 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 434 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 435 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 436 compress = PETSC_FALSE; 437 for (p = pStart; p < pEnd; p++) { 438 PetscInt dof, off, count, offNew, dofNew; 439 440 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 441 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 442 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 443 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 444 count = 0; 445 for (i = 0; i < dof; i++) { 446 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 447 448 if (q >= pStart && q < pEnd) { 449 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 450 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 451 } 452 if (qDof) { 453 PetscInt oldCount = count; 454 455 for (j = 0; j < qDof; j++) { 456 PetscInt k, r = vals[qOff + j]; 457 458 for (k = 0; k < oldCount; k++) { 459 if (valsNew[offNew + k] == r) { 460 break; 461 } 462 } 463 if (k == oldCount) { 464 valsNew[offNew + count++] = r; 465 } 466 } 467 } 468 else { 469 PetscInt k, oldCount = count; 470 471 for (k = 0; k < oldCount; k++) { 472 if (valsNew[offNew + k] == q) { 473 break; 474 } 475 } 476 if (k == oldCount) { 477 valsNew[offNew + count++] = q; 478 } 479 } 480 } 481 if (count < dofNew) { 482 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 483 compress = PETSC_TRUE; 484 } 485 } 486 } 487 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 488 ierr = MPI_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 489 if (!globalAnyNew) { 490 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 491 *sectionNew = NULL; 492 *isNew = NULL; 493 } 494 else { 495 PetscBool globalCompress; 496 497 ierr = MPI_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 498 if (compress) { 499 PetscSection secComp; 500 PetscInt *valsComp = NULL; 501 502 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 503 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 504 for (p = pStart; p < pEnd; p++) { 505 PetscInt dof; 506 507 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 508 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 509 } 510 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 511 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 512 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 513 for (p = pStart; p < pEnd; p++) { 514 PetscInt dof, off, offNew, j; 515 516 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 517 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 518 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 519 for (j = 0; j < dof; j++) { 520 valsComp[offNew + j] = valsNew[off + j]; 521 } 522 } 523 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 524 secNew = secComp; 525 ierr = PetscFree(valsNew);CHKERRQ(ierr); 526 valsNew = valsComp; 527 } 528 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 529 } 530 PetscFunctionReturn(0); 531 } 532 533 #undef __FUNCT__ 534 #define __FUNCT__ "DMPlexComputeConstraints_Tree" 535 static PetscErrorCode DMPlexComputeConstraints_Tree(DM dm) 536 { 537 PetscInt p, pStart, pEnd, *anchors, size; 538 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 539 PetscSection aSec; 540 DMLabel canonLabel; 541 IS aIS; 542 PetscErrorCode ierr; 543 544 PetscFunctionBegin; 545 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 546 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 547 ierr = DMPlexGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr); 548 for (p = pStart; p < pEnd; p++) { 549 PetscInt parent; 550 551 if (canonLabel) { 552 PetscInt canon; 553 554 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 555 if (p != canon) continue; 556 } 557 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 558 if (parent != p) { 559 aMin = PetscMin(aMin,p); 560 aMax = PetscMax(aMax,p+1); 561 } 562 } 563 if (aMin > aMax) { 564 aMin = -1; 565 aMax = -1; 566 } 567 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&aSec);CHKERRQ(ierr); 568 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 569 for (p = aMin; p < aMax; p++) { 570 PetscInt parent, ancestor = p; 571 572 if (canonLabel) { 573 PetscInt canon; 574 575 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 576 if (p != canon) continue; 577 } 578 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 579 while (parent != ancestor) { 580 ancestor = parent; 581 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 582 } 583 if (ancestor != p) { 584 PetscInt closureSize, *closure = NULL; 585 586 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 587 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 588 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 589 } 590 } 591 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 592 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 593 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 594 for (p = aMin; p < aMax; p++) { 595 PetscInt parent, ancestor = p; 596 597 if (canonLabel) { 598 PetscInt canon; 599 600 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 601 if (p != canon) continue; 602 } 603 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 604 while (parent != ancestor) { 605 ancestor = parent; 606 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 607 } 608 if (ancestor != p) { 609 PetscInt j, closureSize, *closure = NULL, aOff; 610 611 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 612 613 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 614 for (j = 0; j < closureSize; j++) { 615 anchors[aOff + j] = closure[2*j]; 616 } 617 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 618 } 619 } 620 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 621 { 622 PetscSection aSecNew = aSec; 623 IS aISNew = aIS; 624 625 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 626 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 627 while (aSecNew) { 628 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 629 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 630 aSec = aSecNew; 631 aIS = aISNew; 632 aSecNew = NULL; 633 aISNew = NULL; 634 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 635 } 636 } 637 ierr = DMPlexSetConstraints(dm,aSec,aIS);CHKERRQ(ierr); 638 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 639 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "DMPlexSetTree_Internal" 645 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical) 646 { 647 DM_Plex *mesh = (DM_Plex *)dm->data; 648 DM refTree; 649 PetscInt size; 650 PetscErrorCode ierr; 651 652 PetscFunctionBegin; 653 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 654 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 655 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 656 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 657 mesh->parentSection = parentSection; 658 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 659 if (parents != mesh->parents) { 660 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 661 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 662 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 663 } 664 if (childIDs != mesh->childIDs) { 665 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 666 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 667 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 668 } 669 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 670 if (refTree) { 671 DMLabel canonLabel; 672 673 ierr = DMPlexGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr); 674 if (canonLabel) { 675 PetscInt i; 676 677 for (i = 0; i < size; i++) { 678 PetscInt canon; 679 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr); 680 if (canon >= 0) { 681 mesh->childIDs[i] = canon; 682 } 683 } 684 } 685 } 686 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 687 if (computeCanonical) { 688 PetscInt d, dim; 689 690 /* add the canonical label */ 691 ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr); 692 ierr = DMPlexCreateLabel(dm,"canonical");CHKERRQ(ierr); 693 for (d = 0; d <= dim; d++) { 694 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 695 const PetscInt *cChildren; 696 697 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 698 for (p = dStart; p < dEnd; p++) { 699 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr); 700 if (cNumChildren) { 701 canon = p; 702 break; 703 } 704 } 705 if (canon == -1) continue; 706 for (p = dStart; p < dEnd; p++) { 707 PetscInt numChildren, i; 708 const PetscInt *children; 709 710 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 711 if (numChildren) { 712 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); 713 ierr = DMPlexSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr); 714 for (i = 0; i < numChildren; i++) { 715 ierr = DMPlexSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr); 716 } 717 } 718 } 719 } 720 } 721 ierr = DMPlexComputeConstraints_Tree(dm);CHKERRQ(ierr); 722 PetscFunctionReturn(0); 723 } 724 725 #undef __FUNCT__ 726 #define __FUNCT__ "DMPlexSetTree" 727 /*@ 728 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 729 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 730 tree root. 731 732 Collective on dm 733 734 Input Parameters: 735 + dm - the DMPlex object 736 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 737 offset indexes the parent and childID list; the reference count of parentSection is incremented 738 . parents - a list of the point parents; copied, can be destroyed 739 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 740 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 741 742 Level: intermediate 743 744 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 745 @*/ 746 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 747 { 748 PetscErrorCode ierr; 749 750 PetscFunctionBegin; 751 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE);CHKERRQ(ierr); 752 PetscFunctionReturn(0); 753 } 754 755 #undef __FUNCT__ 756 #define __FUNCT__ "DMPlexGetTree" 757 /*@ 758 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 759 Collective on dm 760 761 Input Parameters: 762 . dm - the DMPlex object 763 764 Output Parameters: 765 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 766 offset indexes the parent and childID list 767 . parents - a list of the point parents 768 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 769 the child corresponds to the point in the reference tree with index childID 770 . childSection - the inverse of the parent section 771 - children - a list of the point children 772 773 Level: intermediate 774 775 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 776 @*/ 777 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 778 { 779 DM_Plex *mesh = (DM_Plex *)dm->data; 780 781 PetscFunctionBegin; 782 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 783 if (parentSection) *parentSection = mesh->parentSection; 784 if (parents) *parents = mesh->parents; 785 if (childIDs) *childIDs = mesh->childIDs; 786 if (childSection) *childSection = mesh->childSection; 787 if (children) *children = mesh->children; 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNCT__ 792 #define __FUNCT__ "DMPlexGetTreeParent" 793 /*@ 794 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 795 796 Input Parameters: 797 + dm - the DMPlex object 798 - point - the query point 799 800 Output Parameters: 801 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 802 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 803 does not have a parent 804 805 Level: intermediate 806 807 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 808 @*/ 809 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 810 { 811 DM_Plex *mesh = (DM_Plex *)dm->data; 812 PetscSection pSec; 813 PetscErrorCode ierr; 814 815 PetscFunctionBegin; 816 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 817 pSec = mesh->parentSection; 818 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 819 PetscInt dof; 820 821 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 822 if (dof) { 823 PetscInt off; 824 825 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 826 if (parent) *parent = mesh->parents[off]; 827 if (childID) *childID = mesh->childIDs[off]; 828 PetscFunctionReturn(0); 829 } 830 } 831 if (parent) { 832 *parent = point; 833 } 834 if (childID) { 835 *childID = 0; 836 } 837 PetscFunctionReturn(0); 838 } 839 840 #undef __FUNCT__ 841 #define __FUNCT__ "DMPlexGetTreeChildren" 842 /*@C 843 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 844 845 Input Parameters: 846 + dm - the DMPlex object 847 - point - the query point 848 849 Output Parameters: 850 + numChildren - if not NULL, set to the number of children 851 - children - if not NULL, set to a list children, or set to NULL if the point has no children 852 853 Level: intermediate 854 855 Fortran Notes: 856 Since it returns an array, this routine is only available in Fortran 90, and you must 857 include petsc.h90 in your code. 858 859 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 860 @*/ 861 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 862 { 863 DM_Plex *mesh = (DM_Plex *)dm->data; 864 PetscSection childSec; 865 PetscInt dof = 0; 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 870 childSec = mesh->childSection; 871 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 872 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 873 } 874 if (numChildren) *numChildren = dof; 875 if (children) { 876 if (dof) { 877 PetscInt off; 878 879 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 880 *children = &mesh->children[off]; 881 } 882 else { 883 *children = NULL; 884 } 885 } 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "DMPlexComputeConstraintMatrix_ReferenceTree" 891 PetscErrorCode DMPlexComputeConstraintMatrix_ReferenceTree(DM dm) 892 { 893 PetscDS ds; 894 PetscInt spdim; 895 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 896 const PetscInt *anchors; 897 PetscSection section, cSec, aSec; 898 Mat cMat; 899 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 900 IS aIS; 901 PetscErrorCode ierr; 902 903 PetscFunctionBegin; 904 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 905 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 906 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 907 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 908 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 909 ierr = DMPlexGetConstraintSection(dm,&cSec);CHKERRQ(ierr); 910 ierr = DMPlexGetConstraints(dm,&aSec,&aIS);CHKERRQ(ierr); 911 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 912 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 913 ierr = DMPlexGetConstraintMatrix(dm,&cMat);CHKERRQ(ierr); 914 ierr = DMPlexGetDimension(dm,&spdim);CHKERRQ(ierr); 915 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 916 917 for (f = 0; f < numFields; f++) { 918 PetscFE fe; 919 PetscDualSpace space; 920 PetscInt i, j, k, nPoints, offset; 921 PetscInt fSize, fComp; 922 PetscScalar *B = NULL; 923 PetscReal *weights, *pointsRef, *pointsReal; 924 Mat Amat, Bmat, Xmat; 925 926 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *)(&fe));CHKERRQ(ierr); 927 ierr = PetscFEGetDualSpace(fe,&space);CHKERRQ(ierr); 928 ierr = PetscDualSpaceGetDimension(space,&fSize);CHKERRQ(ierr); 929 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 930 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 931 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 932 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 933 ierr = MatSetUp(Amat);CHKERRQ(ierr); 934 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 935 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 936 nPoints = 0; 937 for (i = 0; i < fSize; i++) { 938 PetscInt qPoints; 939 PetscQuadrature quad; 940 941 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 942 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 943 nPoints += qPoints; 944 } 945 ierr = PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);CHKERRQ(ierr); 946 offset = 0; 947 for (i = 0; i < fSize; i++) { 948 PetscInt qPoints; 949 const PetscReal *p, *w; 950 PetscQuadrature quad; 951 952 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 953 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 954 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 955 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 956 offset += qPoints; 957 } 958 ierr = PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 959 offset = 0; 960 for (i = 0; i < fSize; i++) { 961 PetscInt qPoints; 962 PetscQuadrature quad; 963 964 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 965 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 966 for (j = 0; j < fSize; j++) { 967 PetscScalar val = 0.; 968 969 for (k = 0; k < qPoints; k++) { 970 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 971 } 972 ierr = MatSetValue(Amat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 973 } 974 offset += qPoints; 975 } 976 ierr = PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 977 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 978 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 979 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 980 for (c = cStart; c < cEnd; c++) { 981 PetscInt parent; 982 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 983 PetscInt *childOffsets, *parentOffsets; 984 985 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 986 if (parent == c) continue; 987 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 988 for (i = 0; i < closureSize; i++) { 989 PetscInt p = closure[2*i]; 990 PetscInt conDof; 991 992 if (p < conStart || p >= conEnd) continue; 993 if (numFields > 1) { 994 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 995 } 996 else { 997 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 998 } 999 if (conDof) break; 1000 } 1001 if (i == closureSize) { 1002 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1003 continue; 1004 } 1005 1006 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 1007 ierr = DMPlexComputeCellGeometry(dm, parent, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 1008 for (i = 0; i < nPoints; i++) { 1009 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);CHKERRQ(ierr); 1010 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);CHKERRQ(ierr); 1011 } 1012 ierr = PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1013 offset = 0; 1014 for (i = 0; i < fSize; i++) { 1015 PetscInt qPoints; 1016 PetscQuadrature quad; 1017 1018 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1019 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1020 for (j = 0; j < fSize; j++) { 1021 PetscScalar val = 0.; 1022 1023 for (k = 0; k < qPoints; k++) { 1024 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 1025 } 1026 MatSetValue(Bmat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 1027 } 1028 offset += qPoints; 1029 } 1030 ierr = PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1031 ierr = MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1032 ierr = MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1033 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1034 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1035 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1036 childOffsets[0] = 0; 1037 for (i = 0; i < closureSize; i++) { 1038 PetscInt p = closure[2*i]; 1039 PetscInt dof; 1040 1041 if (numFields > 1) { 1042 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1043 } 1044 else { 1045 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1046 } 1047 childOffsets[i+1]=childOffsets[i]+dof / fComp; 1048 } 1049 parentOffsets[0] = 0; 1050 for (i = 0; i < closureSizeP; i++) { 1051 PetscInt p = closureP[2*i]; 1052 PetscInt dof; 1053 1054 if (numFields > 1) { 1055 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1056 } 1057 else { 1058 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1059 } 1060 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 1061 } 1062 for (i = 0; i < closureSize; i++) { 1063 PetscInt conDof, conOff, aDof, aOff; 1064 PetscInt p = closure[2*i]; 1065 PetscInt o = closure[2*i+1]; 1066 1067 if (p < conStart || p >= conEnd) continue; 1068 if (numFields > 1) { 1069 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1070 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1071 } 1072 else { 1073 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1074 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1075 } 1076 if (!conDof) continue; 1077 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1078 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1079 for (k = 0; k < aDof; k++) { 1080 PetscInt a = anchors[aOff + k]; 1081 PetscInt aSecDof, aSecOff; 1082 1083 if (numFields > 1) { 1084 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1085 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1086 } 1087 else { 1088 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1089 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1090 } 1091 if (!aSecDof) continue; 1092 1093 for (j = 0; j < closureSizeP; j++) { 1094 PetscInt q = closureP[2*j]; 1095 PetscInt oq = closureP[2*j+1]; 1096 1097 if (q == a) { 1098 PetscInt r, s, t; 1099 1100 for (r = childOffsets[i]; r < childOffsets[i+1]; r++) { 1101 for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) { 1102 PetscScalar val; 1103 PetscInt insertCol, insertRow; 1104 1105 ierr = MatGetValue(Xmat,r,s,&val);CHKERRQ(ierr); 1106 if (o >= 0) { 1107 insertRow = conOff + fComp * (r - childOffsets[i]); 1108 } 1109 else { 1110 insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r); 1111 } 1112 if (oq >= 0) { 1113 insertCol = aSecOff + fComp * (s - parentOffsets[j]); 1114 } 1115 else { 1116 insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s); 1117 } 1118 for (t = 0; t < fComp; t++) { 1119 ierr = MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);CHKERRQ(ierr); 1120 } 1121 } 1122 } 1123 } 1124 } 1125 } 1126 } 1127 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1128 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1129 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1130 } 1131 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1132 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1133 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1134 ierr = PetscFree3(weights,pointsRef,pointsReal);CHKERRQ(ierr); 1135 } 1136 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1137 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1138 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1139 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1140 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "DMPlexComputeConstraintMatrix_Tree" 1146 PetscErrorCode DMPlexComputeConstraintMatrix_Tree(DM dm) 1147 { 1148 DM refTree; 1149 PetscDS ds; 1150 Mat refCmat, cMat; 1151 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1152 PetscScalar ***refPointFieldMats, *pointWork; 1153 PetscSection refConSec, conSec, refAnSec, anSec, section, refSection; 1154 IS refAnIS, anIS; 1155 const PetscInt *refAnchors, *anchors; 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1160 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1161 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1162 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 1163 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 1164 ierr = DMPlexComputeConstraintMatrix_ReferenceTree(refTree);CHKERRQ(ierr); 1165 ierr = DMPlexGetConstraintMatrix(refTree,&refCmat);CHKERRQ(ierr); 1166 ierr = DMPlexGetChart(refTree,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1167 ierr = DMPlexGetConstraints(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1168 ierr = DMPlexGetConstraints(dm,&anSec,&anIS);CHKERRQ(ierr); 1169 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1170 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr); 1171 ierr = DMPlexGetConstraintSection(refTree,&refConSec);CHKERRQ(ierr); 1172 ierr = DMPlexGetConstraintSection(dm,&conSec);CHKERRQ(ierr); 1173 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 1174 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 1175 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1176 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr); 1177 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 1178 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr); 1179 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1180 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1181 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 1182 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr); 1183 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr); 1184 1185 /* step 1: get submats for every constrained point in the reference tree */ 1186 for (p = pRefStart; p < pRefEnd; p++) { 1187 PetscInt parent, closureSize, *closure = NULL, pDof; 1188 1189 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1190 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1191 if (!pDof || parent == p) continue; 1192 1193 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 1194 ierr = PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr); 1195 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1196 for (f = 0; f < numFields; f++) { 1197 PetscInt cDof, cOff, numCols, r, i, fComp; 1198 PetscFE fe; 1199 1200 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);CHKERRQ(ierr); 1201 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1202 1203 if (numFields > 1) { 1204 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1205 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 1206 } 1207 else { 1208 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1209 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 1210 } 1211 1212 if (!cDof) continue; 1213 for (r = 0; r < cDof; r++) { 1214 rows[r] = cOff + r; 1215 } 1216 numCols = 0; 1217 for (i = 0; i < closureSize; i++) { 1218 PetscInt q = closure[2*i]; 1219 PetscInt o = closure[2*i+1]; 1220 PetscInt aDof, aOff, j; 1221 1222 if (numFields > 1) { 1223 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1224 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1225 } 1226 else { 1227 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1228 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1229 } 1230 1231 for (j = 0; j < aDof; j++) { 1232 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1233 PetscInt comp = (j % fComp); 1234 1235 cols[numCols++] = aOff + node * fComp + comp; 1236 } 1237 } 1238 refPointFieldN[p-pRefStart][f] = numCols; 1239 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1240 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1241 } 1242 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1243 } 1244 1245 /* step 2: compute the preorder */ 1246 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1247 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr); 1248 for (p = pStart; p < pEnd; p++) { 1249 perm[p - pStart] = p; 1250 iperm[p - pStart] = p-pStart; 1251 } 1252 for (p = 0; p < pEnd - pStart;) { 1253 PetscInt point = perm[p]; 1254 PetscInt parent; 1255 1256 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr); 1257 if (parent == point) { 1258 p++; 1259 } 1260 else { 1261 PetscInt size, closureSize, *closure = NULL, i; 1262 1263 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1264 for (i = 0; i < closureSize; i++) { 1265 PetscInt q = closure[2*i]; 1266 if (iperm[q-pStart] > iperm[point-pStart]) { 1267 /* swap */ 1268 perm[p] = q; 1269 perm[iperm[q-pStart]] = point; 1270 iperm[point-pStart] = iperm[q-pStart]; 1271 iperm[q-pStart] = p; 1272 break; 1273 } 1274 } 1275 size = closureSize; 1276 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1277 if (i == size) { 1278 p++; 1279 } 1280 } 1281 } 1282 1283 /* step 3: fill the constraint matrix */ 1284 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1285 * allow progressive fill without assembly, so we are going to set up the 1286 * values outside of the Mat first. 1287 */ 1288 { 1289 PetscInt nRows, row, nnz; 1290 PetscBool done; 1291 const PetscInt *ia, *ja; 1292 PetscScalar *vals; 1293 1294 ierr = DMPlexGetConstraintMatrix(dm,&cMat);CHKERRQ(ierr); 1295 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1296 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1297 nnz = ia[nRows]; 1298 /* malloc and then zero rows right before we fill them: this way valgrind 1299 * can tell if we are doing progressive fill in the wrong order */ 1300 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr); 1301 for (p = 0; p < pEnd - pStart; p++) { 1302 PetscInt parent, childid, closureSize, *closure = NULL; 1303 PetscInt point = perm[p], pointDof; 1304 1305 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr); 1306 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1307 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr); 1308 if (!pointDof) continue; 1309 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1310 for (f = 0; f < numFields; f++) { 1311 PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp, matOffset, offset; 1312 PetscScalar *pointMat; 1313 PetscFE fe; 1314 1315 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);CHKERRQ(ierr); 1316 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1317 1318 if (numFields > 1) { 1319 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr); 1320 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr); 1321 } 1322 else { 1323 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr); 1324 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr); 1325 } 1326 if (!cDof) continue; 1327 1328 /* make sure that every row for this point is the same size */ 1329 #if defined(PETSC_USE_DEBUG) 1330 for (r = 0; r < cDof; r++) { 1331 if (cDof > 1 && r) { 1332 if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) { 1333 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %d vs. %d", (ia[rows[r]+1]-ia[rows[r]]), (ia[rows[r]]-ia[rows[r]-1])); 1334 } 1335 } 1336 } 1337 #endif 1338 /* zero rows */ 1339 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1340 vals[i] = 0.; 1341 } 1342 matOffset = ia[cOff]; 1343 numFillCols = ia[cOff+1] - matOffset; 1344 pointMat = refPointFieldMats[childid-pRefStart][f]; 1345 numCols = refPointFieldN[childid-pRefStart][f]; 1346 offset = 0; 1347 for (i = 0; i < closureSize; i++) { 1348 PetscInt q = closure[2*i]; 1349 PetscInt o = closure[2*i+1]; 1350 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1351 1352 qConDof = qConOff = 0; 1353 if (numFields > 1) { 1354 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr); 1355 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr); 1356 if (q >= conStart && q < conEnd) { 1357 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr); 1358 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr); 1359 } 1360 } 1361 else { 1362 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr); 1363 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr); 1364 if (q >= conStart && q < conEnd) { 1365 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr); 1366 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr); 1367 } 1368 } 1369 if (!aDof) continue; 1370 if (qConDof) { 1371 /* this point has anchors: its rows of the matrix should already 1372 * be filled, thanks to preordering */ 1373 /* first multiply into pointWork, then set in matrix */ 1374 PetscInt aMatOffset = ia[qConOff]; 1375 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1376 for (r = 0; r < cDof; r++) { 1377 for (j = 0; j < aNumFillCols; j++) { 1378 PetscScalar inVal = 0; 1379 for (k = 0; k < aDof; k++) { 1380 PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp); 1381 PetscInt comp = (k % fComp); 1382 PetscInt col = node * fComp + comp; 1383 1384 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j]; 1385 } 1386 pointWork[r * aNumFillCols + j] = inVal; 1387 } 1388 } 1389 /* assume that the columns are sorted, spend less time searching */ 1390 for (j = 0, k = 0; j < aNumFillCols; j++) { 1391 PetscInt col = ja[aMatOffset + j]; 1392 for (;k < numFillCols; k++) { 1393 if (ja[matOffset + k] == col) { 1394 break; 1395 } 1396 } 1397 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col); 1398 for (r = 0; r < cDof; r++) { 1399 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1400 } 1401 } 1402 } 1403 else { 1404 /* find where to put this portion of pointMat into the matrix */ 1405 for (k = 0; k < numFillCols; k++) { 1406 if (ja[matOffset + k] == aOff) { 1407 break; 1408 } 1409 } 1410 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff); 1411 for (j = 0; j < aDof; j++) { 1412 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1413 PetscInt comp = (j % fComp); 1414 PetscInt col = node * fComp + comp; 1415 for (r = 0; r < cDof; r++) { 1416 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col]; 1417 } 1418 } 1419 } 1420 offset += aDof; 1421 } 1422 } 1423 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1424 } 1425 for (row = 0; row < nRows; row++) { 1426 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr); 1427 } 1428 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1429 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1430 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1431 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1432 ierr = PetscFree(vals);CHKERRQ(ierr); 1433 } 1434 1435 /* clean up */ 1436 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1437 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr); 1438 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr); 1439 ierr = PetscFree(rows);CHKERRQ(ierr); 1440 ierr = PetscFree(cols);CHKERRQ(ierr); 1441 ierr = PetscFree(pointWork);CHKERRQ(ierr); 1442 for (p = pRefStart; p < pRefEnd; p++) { 1443 PetscInt parent, pDof; 1444 1445 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1446 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1447 if (!pDof || parent == p) continue; 1448 1449 for (f = 0; f < numFields; f++) { 1450 PetscInt cDof; 1451 1452 if (numFields > 1) { 1453 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1454 } 1455 else { 1456 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1457 } 1458 1459 if (!cDof) continue; 1460 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 1461 } 1462 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 1463 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr); 1464 } 1465 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 1466 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr); 1467 PetscFunctionReturn(0); 1468 } 1469 1470