1 #include <petsc/private/pcpatchimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* For ksp->setfromoptionscalled */ 3 #include <petsc/private/dmpleximpl.h> /* For DMPlexComputeJacobian_Patch_Internal() */ 4 #include <petscsf.h> 5 #include <petscbt.h> 6 #include <petscds.h> 7 8 PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Scatter, PC_Patch_Apply, PC_Patch_Prealloc; 9 10 PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format) 11 { 12 PetscErrorCode ierr; 13 14 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 15 ierr = PetscObjectView(obj, viewer);CHKERRQ(ierr); 16 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 17 return(0); 18 } 19 20 static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 21 { 22 PetscInt starSize; 23 PetscInt *star = NULL, si; 24 PetscErrorCode ierr; 25 26 PetscFunctionBegin; 27 PetscHSetIClear(ht); 28 /* To start with, add the point we care about */ 29 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 30 /* Loop over all the points that this point connects to */ 31 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 32 for (si = 0; si < starSize*2; si += 2) {ierr = PetscHSetIAdd(ht, star[si]);CHKERRQ(ierr);} 33 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 38 { 39 PC_PATCH *patch = (PC_PATCH *) vpatch; 40 PetscInt starSize; 41 PetscInt *star = NULL; 42 PetscBool shouldIgnore = PETSC_FALSE; 43 PetscInt cStart, cEnd, iStart, iEnd, si; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 48 /* To start with, add the point we care about */ 49 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 50 /* Should we ignore any points of a certain dimension? */ 51 if (patch->vankadim >= 0) { 52 shouldIgnore = PETSC_TRUE; 53 ierr = DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd);CHKERRQ(ierr); 54 } 55 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 56 /* Loop over all the cells that this point connects to */ 57 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 58 for (si = 0; si < starSize*2; si += 2) { 59 const PetscInt cell = star[si]; 60 PetscInt closureSize; 61 PetscInt *closure = NULL, ci; 62 63 if (cell < cStart || cell >= cEnd) continue; 64 /* now loop over all entities in the closure of that cell */ 65 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 66 for (ci = 0; ci < closureSize*2; ci += 2) { 67 const PetscInt newpoint = closure[ci]; 68 69 /* We've been told to ignore entities of this type.*/ 70 if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue; 71 ierr = PetscHSetIAdd(ht, newpoint);CHKERRQ(ierr); 72 } 73 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 74 } 75 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 80 { 81 PC_PATCH *patch = (PC_PATCH *) vpatch; 82 DMLabel ghost = NULL; 83 const PetscInt *leaves; 84 PetscInt nleaves, pStart, pEnd, loc; 85 PetscBool isFiredrake; 86 DM plex; 87 PetscBool flg; 88 PetscInt starSize; 89 PetscInt *star = NULL; 90 PetscInt opoint, overlapi; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 PetscHSetIClear(ht); 95 96 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 97 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 98 99 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 100 if (isFiredrake) { 101 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 102 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 103 } else { 104 PetscSF sf; 105 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 106 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 107 nleaves = PetscMax(nleaves, 0); 108 } 109 110 for (opoint = pStart; opoint < pEnd; ++opoint) { 111 if (ghost) {ierr = DMLabelHasPoint(ghost, opoint, &flg);CHKERRQ(ierr);} 112 else {ierr = PetscFindInt(opoint, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 113 /* Not an owned entity, don't make a cell patch. */ 114 if (flg) continue; 115 ierr = PetscHSetIAdd(ht, opoint);CHKERRQ(ierr); 116 } 117 118 /* Now build the overlap for the patch */ 119 for (overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) { 120 PetscInt index = 0; 121 PetscInt *htpoints = NULL; 122 PetscInt htsize; 123 PetscInt i; 124 125 ierr = PetscHSetIGetSize(ht, &htsize);CHKERRQ(ierr); 126 ierr = PetscMalloc1(htsize, &htpoints);CHKERRQ(ierr); 127 ierr = PetscHSetIGetElems(ht, &index, htpoints);CHKERRQ(ierr); 128 129 for (i = 0; i < htsize; ++i) { 130 PetscInt hpoint = htpoints[i]; 131 PetscInt si; 132 133 ierr = DMPlexGetTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 134 for (si = 0; si < starSize*2; si += 2) { 135 const PetscInt starp = star[si]; 136 PetscInt closureSize; 137 PetscInt *closure = NULL, ci; 138 139 /* now loop over all entities in the closure of starp */ 140 ierr = DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 141 for (ci = 0; ci < closureSize*2; ci += 2) { 142 const PetscInt closstarp = closure[ci]; 143 ierr = PetscHSetIAdd(ht, closstarp);CHKERRQ(ierr); 144 } 145 ierr = DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 146 } 147 ierr = DMPlexRestoreTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 148 } 149 ierr = PetscFree(htpoints);CHKERRQ(ierr); 150 } 151 152 PetscFunctionReturn(0); 153 } 154 155 /* The user's already set the patches in patch->userIS. Build the hash tables */ 156 static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 157 { 158 PC_PATCH *patch = (PC_PATCH *) vpatch; 159 IS patchis = patch->userIS[point]; 160 PetscInt n; 161 const PetscInt *patchdata; 162 PetscInt pStart, pEnd, i; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 167 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 168 ierr = ISGetLocalSize(patchis, &n);CHKERRQ(ierr); 169 ierr = ISGetIndices(patchis, &patchdata);CHKERRQ(ierr); 170 for (i = 0; i < n; ++i) { 171 const PetscInt ownedpoint = patchdata[i]; 172 173 if (ownedpoint < pStart || ownedpoint >= pEnd) { 174 SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd); 175 } 176 ierr = PetscHSetIAdd(ht, ownedpoint);CHKERRQ(ierr); 177 } 178 ierr = ISRestoreIndices(patchis, &patchdata);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs) 183 { 184 PC_PATCH *patch = (PC_PATCH *) pc->data; 185 PetscInt i; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 if (n == 1 && bs[0] == 1) { 190 patch->defaultSF = sf[0]; 191 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 192 } else { 193 PetscInt allRoots = 0, allLeaves = 0; 194 PetscInt leafOffset = 0; 195 PetscInt *ilocal = NULL; 196 PetscSFNode *iremote = NULL; 197 PetscInt *remoteOffsets = NULL; 198 PetscInt index = 0; 199 PetscHMapI rankToIndex; 200 PetscInt numRanks = 0; 201 PetscSFNode *remote = NULL; 202 PetscSF rankSF; 203 PetscInt *ranks = NULL; 204 PetscInt *offsets = NULL; 205 MPI_Datatype contig; 206 PetscHSetI ranksUniq; 207 208 /* First figure out how many dofs there are in the concatenated numbering. 209 * allRoots: number of owned global dofs; 210 * allLeaves: number of visible dofs (global + ghosted). 211 */ 212 for (i = 0; i < n; ++i) { 213 PetscInt nroots, nleaves; 214 215 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr); 216 allRoots += nroots * bs[i]; 217 allLeaves += nleaves * bs[i]; 218 } 219 ierr = PetscMalloc1(allLeaves, &ilocal);CHKERRQ(ierr); 220 ierr = PetscMalloc1(allLeaves, &iremote);CHKERRQ(ierr); 221 /* Now build an SF that just contains process connectivity. */ 222 ierr = PetscHSetICreate(&ranksUniq);CHKERRQ(ierr); 223 for (i = 0; i < n; ++i) { 224 const PetscMPIInt *ranks = NULL; 225 PetscInt nranks, j; 226 227 ierr = PetscSFSetUp(sf[i]);CHKERRQ(ierr); 228 ierr = PetscSFGetRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL);CHKERRQ(ierr); 229 /* These are all the ranks who communicate with me. */ 230 for (j = 0; j < nranks; ++j) { 231 ierr = PetscHSetIAdd(ranksUniq, (PetscInt) ranks[j]);CHKERRQ(ierr); 232 } 233 } 234 ierr = PetscHSetIGetSize(ranksUniq, &numRanks);CHKERRQ(ierr); 235 ierr = PetscMalloc1(numRanks, &remote);CHKERRQ(ierr); 236 ierr = PetscMalloc1(numRanks, &ranks);CHKERRQ(ierr); 237 ierr = PetscHSetIGetElems(ranksUniq, &index, ranks);CHKERRQ(ierr); 238 239 ierr = PetscHMapICreate(&rankToIndex);CHKERRQ(ierr); 240 for (i = 0; i < numRanks; ++i) { 241 remote[i].rank = ranks[i]; 242 remote[i].index = 0; 243 ierr = PetscHMapISet(rankToIndex, ranks[i], i);CHKERRQ(ierr); 244 } 245 ierr = PetscFree(ranks);CHKERRQ(ierr); 246 ierr = PetscHSetIDestroy(&ranksUniq);CHKERRQ(ierr); 247 ierr = PetscSFCreate(PetscObjectComm((PetscObject) pc), &rankSF);CHKERRQ(ierr); 248 ierr = PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 249 ierr = PetscSFSetUp(rankSF);CHKERRQ(ierr); 250 /* OK, use it to communicate the root offset on the remote 251 * processes for each subspace. */ 252 ierr = PetscMalloc1(n, &offsets);CHKERRQ(ierr); 253 ierr = PetscMalloc1(n*numRanks, &remoteOffsets);CHKERRQ(ierr); 254 255 offsets[0] = 0; 256 for (i = 1; i < n; ++i) { 257 PetscInt nroots; 258 259 ierr = PetscSFGetGraph(sf[i-1], &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 260 offsets[i] = offsets[i-1] + nroots*bs[i-1]; 261 } 262 /* Offsets are the offsets on the current process of the 263 * global dof numbering for the subspaces. */ 264 ierr = MPI_Type_contiguous(n, MPIU_INT, &contig);CHKERRQ(ierr); 265 ierr = MPI_Type_commit(&contig);CHKERRQ(ierr); 266 267 ierr = PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 268 ierr = PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 269 ierr = MPI_Type_free(&contig);CHKERRQ(ierr); 270 ierr = PetscFree(offsets);CHKERRQ(ierr); 271 ierr = PetscSFDestroy(&rankSF);CHKERRQ(ierr); 272 /* Now remoteOffsets contains the offsets on the remote 273 * processes who communicate with me. So now we can 274 * concatenate the list of SFs into a single one. */ 275 index = 0; 276 for (i = 0; i < n; ++i) { 277 const PetscSFNode *remote = NULL; 278 const PetscInt *local = NULL; 279 PetscInt nroots, nleaves, j; 280 281 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);CHKERRQ(ierr); 282 for (j = 0; j < nleaves; ++j) { 283 PetscInt rank = remote[j].rank; 284 PetscInt idx, rootOffset, k; 285 286 ierr = PetscHMapIGet(rankToIndex, rank, &idx);CHKERRQ(ierr); 287 if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?"); 288 /* Offset on given rank for ith subspace */ 289 rootOffset = remoteOffsets[n*idx + i]; 290 for (k = 0; k < bs[i]; ++k) { 291 ilocal[index] = (local ? local[j] : j)*bs[i] + k + leafOffset; 292 iremote[index].rank = remote[j].rank; 293 iremote[index].index = remote[j].index*bs[i] + k + rootOffset; 294 ++index; 295 } 296 } 297 leafOffset += nleaves * bs[i]; 298 } 299 ierr = PetscHMapIDestroy(&rankToIndex);CHKERRQ(ierr); 300 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 301 ierr = PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->defaultSF);CHKERRQ(ierr); 302 ierr = PetscSFSetGraph(patch->defaultSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);CHKERRQ(ierr); 303 } 304 PetscFunctionReturn(0); 305 } 306 307 /* TODO: Docs */ 308 PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim) 309 { 310 PC_PATCH *patch = (PC_PATCH *) pc->data; 311 PetscFunctionBegin; 312 patch->ignoredim = dim; 313 PetscFunctionReturn(0); 314 } 315 316 /* TODO: Docs */ 317 PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim) 318 { 319 PC_PATCH *patch = (PC_PATCH *) pc->data; 320 PetscFunctionBegin; 321 *dim = patch->ignoredim; 322 PetscFunctionReturn(0); 323 } 324 325 /* TODO: Docs */ 326 PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg) 327 { 328 PC_PATCH *patch = (PC_PATCH *) pc->data; 329 PetscFunctionBegin; 330 patch->save_operators = flg; 331 PetscFunctionReturn(0); 332 } 333 334 /* TODO: Docs */ 335 PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg) 336 { 337 PC_PATCH *patch = (PC_PATCH *) pc->data; 338 PetscFunctionBegin; 339 *flg = patch->save_operators; 340 PetscFunctionReturn(0); 341 } 342 343 /* TODO: Docs */ 344 PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg) 345 { 346 PC_PATCH *patch = (PC_PATCH *) pc->data; 347 PetscFunctionBegin; 348 patch->precomputeElementTensors = flg; 349 PetscFunctionReturn(0); 350 } 351 352 /* TODO: Docs */ 353 PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg) 354 { 355 PC_PATCH *patch = (PC_PATCH *) pc->data; 356 PetscFunctionBegin; 357 *flg = patch->precomputeElementTensors; 358 PetscFunctionReturn(0); 359 } 360 361 /* TODO: Docs */ 362 PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg) 363 { 364 PC_PATCH *patch = (PC_PATCH *) pc->data; 365 PetscFunctionBegin; 366 patch->partition_of_unity = flg; 367 PetscFunctionReturn(0); 368 } 369 370 /* TODO: Docs */ 371 PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg) 372 { 373 PC_PATCH *patch = (PC_PATCH *) pc->data; 374 PetscFunctionBegin; 375 *flg = patch->partition_of_unity; 376 PetscFunctionReturn(0); 377 } 378 379 /* TODO: Docs */ 380 PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type) 381 { 382 PC_PATCH *patch = (PC_PATCH *) pc->data; 383 PetscFunctionBegin; 384 if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type"); 385 patch->local_composition_type = type; 386 PetscFunctionReturn(0); 387 } 388 389 /* TODO: Docs */ 390 PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type) 391 { 392 PC_PATCH *patch = (PC_PATCH *) pc->data; 393 PetscFunctionBegin; 394 *type = patch->local_composition_type; 395 PetscFunctionReturn(0); 396 } 397 398 /* TODO: Docs */ 399 PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type) 400 { 401 PC_PATCH *patch = (PC_PATCH *) pc->data; 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 if (patch->sub_mat_type) {ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr);} 406 ierr = PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 407 PetscFunctionReturn(0); 408 } 409 410 /* TODO: Docs */ 411 PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type) 412 { 413 PC_PATCH *patch = (PC_PATCH *) pc->data; 414 PetscFunctionBegin; 415 *sub_mat_type = patch->sub_mat_type; 416 PetscFunctionReturn(0); 417 } 418 419 /* TODO: Docs */ 420 PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering) 421 { 422 PC_PATCH *patch = (PC_PATCH *) pc->data; 423 PetscErrorCode ierr; 424 425 PetscFunctionBegin; 426 patch->cellNumbering = cellNumbering; 427 ierr = PetscObjectReference((PetscObject) cellNumbering);CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 431 /* TODO: Docs */ 432 PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering) 433 { 434 PC_PATCH *patch = (PC_PATCH *) pc->data; 435 PetscFunctionBegin; 436 *cellNumbering = patch->cellNumbering; 437 PetscFunctionReturn(0); 438 } 439 440 /* TODO: Docs */ 441 PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 442 { 443 PC_PATCH *patch = (PC_PATCH *) pc->data; 444 445 PetscFunctionBegin; 446 patch->ctype = ctype; 447 switch (ctype) { 448 case PC_PATCH_STAR: 449 patch->user_patches = PETSC_FALSE; 450 patch->patchconstructop = PCPatchConstruct_Star; 451 break; 452 case PC_PATCH_VANKA: 453 patch->user_patches = PETSC_FALSE; 454 patch->patchconstructop = PCPatchConstruct_Vanka; 455 break; 456 case PC_PATCH_PARDECOMP: 457 patch->user_patches = PETSC_FALSE; 458 patch->patchconstructop = PCPatchConstruct_Pardecomp; 459 break; 460 case PC_PATCH_USER: 461 case PC_PATCH_PYTHON: 462 patch->user_patches = PETSC_TRUE; 463 patch->patchconstructop = PCPatchConstruct_User; 464 if (func) { 465 patch->userpatchconstructionop = func; 466 patch->userpatchconstructctx = ctx; 467 } 468 break; 469 default: 470 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 471 } 472 PetscFunctionReturn(0); 473 } 474 475 /* TODO: Docs */ 476 PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx) 477 { 478 PC_PATCH *patch = (PC_PATCH *) pc->data; 479 480 PetscFunctionBegin; 481 *ctype = patch->ctype; 482 switch (patch->ctype) { 483 case PC_PATCH_STAR: 484 case PC_PATCH_VANKA: 485 case PC_PATCH_PARDECOMP: 486 break; 487 case PC_PATCH_USER: 488 case PC_PATCH_PYTHON: 489 *func = patch->userpatchconstructionop; 490 *ctx = patch->userpatchconstructctx; 491 break; 492 default: 493 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 494 } 495 PetscFunctionReturn(0); 496 } 497 498 /* TODO: Docs */ 499 PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 500 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 501 { 502 PC_PATCH *patch = (PC_PATCH *) pc->data; 503 DM dm; 504 PetscSF *sfs; 505 PetscInt cStart, cEnd, i, j; 506 PetscErrorCode ierr; 507 508 PetscFunctionBegin; 509 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 510 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 511 ierr = PetscMalloc1(nsubspaces, &sfs);CHKERRQ(ierr); 512 ierr = PetscMalloc1(nsubspaces, &patch->dofSection);CHKERRQ(ierr); 513 ierr = PetscMalloc1(nsubspaces, &patch->bs);CHKERRQ(ierr); 514 ierr = PetscMalloc1(nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 515 ierr = PetscMalloc1(nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 516 ierr = PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 517 518 patch->nsubspaces = nsubspaces; 519 patch->totalDofsPerCell = 0; 520 for (i = 0; i < nsubspaces; ++i) { 521 ierr = DMGetDefaultSection(dms[i], &patch->dofSection[i]);CHKERRQ(ierr); 522 ierr = PetscObjectReference((PetscObject) patch->dofSection[i]);CHKERRQ(ierr); 523 ierr = DMGetDefaultSF(dms[i], &sfs[i]);CHKERRQ(ierr); 524 patch->bs[i] = bs[i]; 525 patch->nodesPerCell[i] = nodesPerCell[i]; 526 patch->totalDofsPerCell += nodesPerCell[i]*bs[i]; 527 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 528 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 529 patch->subspaceOffsets[i] = subspaceOffsets[i]; 530 } 531 ierr = PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);CHKERRQ(ierr); 532 ierr = PetscFree(sfs);CHKERRQ(ierr); 533 534 patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces]; 535 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 536 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 /* TODO: Docs */ 541 PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 542 { 543 PC_PATCH *patch = (PC_PATCH *) pc->data; 544 PetscInt cStart, cEnd, i, j; 545 PetscErrorCode ierr; 546 547 PetscFunctionBegin; 548 patch->combined = PETSC_TRUE; 549 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 550 ierr = DMGetNumFields(dm, &patch->nsubspaces);CHKERRQ(ierr); 551 ierr = PetscCalloc1(patch->nsubspaces, &patch->dofSection);CHKERRQ(ierr); 552 ierr = PetscMalloc1(patch->nsubspaces, &patch->bs);CHKERRQ(ierr); 553 ierr = PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 554 ierr = PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 555 ierr = PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 556 ierr = DMGetDefaultSection(dm, &patch->dofSection[0]);CHKERRQ(ierr); 557 ierr = PetscObjectReference((PetscObject) patch->dofSection[0]);CHKERRQ(ierr); 558 ierr = PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);CHKERRQ(ierr); 559 patch->totalDofsPerCell = 0; 560 for (i = 0; i < patch->nsubspaces; ++i) { 561 patch->bs[i] = 1; 562 patch->nodesPerCell[i] = nodesPerCell[i]; 563 patch->totalDofsPerCell += nodesPerCell[i]; 564 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 565 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 566 } 567 ierr = DMGetDefaultSF(dm, &patch->defaultSF);CHKERRQ(ierr); 568 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 569 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 570 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 571 PetscFunctionReturn(0); 572 } 573 574 /*@C 575 576 PCPatchSetComputeFunction - Set the callback used to compute patch residuals 577 578 Logically collective on PC 579 580 Input Parameters: 581 + pc - The PC 582 . func - The callback 583 - ctx - The user context 584 585 Calling sequence of func: 586 $ func (PC pc,PetscInt point,Vec x,Vec f,IS cellIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx) 587 588 + pc - The PC 589 . point - The point 590 . x - The input solution (not used in linear problems) 591 . f - The patch residual vector 592 . cellIS - An array of the cell numbers 593 . n - The size of dofsArray 594 . dofsArray - The dofmap for the dofs to be solved for 595 . dofsArrayWithAll - The dofmap for all dofs on the patch 596 - ctx - The user context 597 598 Level: advanced 599 600 Notes: 601 The matrix entries have been set to zero before the call. 602 603 .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 604 @*/ 605 PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 606 { 607 PC_PATCH *patch = (PC_PATCH *) pc->data; 608 609 PetscFunctionBegin; 610 patch->usercomputef = func; 611 patch->usercomputefctx = ctx; 612 PetscFunctionReturn(0); 613 } 614 615 /*@C 616 617 PCPatchSetComputeFunctionInteriorFacets - Set the callback used to compute facet integrals for patch residuals 618 619 Logically collective on PC 620 621 Input Parameters: 622 + pc - The PC 623 . func - The callback 624 - ctx - The user context 625 626 Calling sequence of func: 627 $ func (PC pc,PetscInt point,Vec x,Vec f,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx) 628 629 + pc - The PC 630 . point - The point 631 . x - The input solution (not used in linear problems) 632 . f - The patch residual vector 633 . facetIS - An array of the facet numbers 634 . n - The size of dofsArray 635 . dofsArray - The dofmap for the dofs to be solved for 636 . dofsArrayWithAll - The dofmap for all dofs on the patch 637 - ctx - The user context 638 639 Level: advanced 640 641 Notes: 642 The matrix entries have been set to zero before the call. 643 644 .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 645 @*/ 646 PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 647 { 648 PC_PATCH *patch = (PC_PATCH *) pc->data; 649 650 PetscFunctionBegin; 651 patch->usercomputefintfacet = func; 652 patch->usercomputefintfacetctx = ctx; 653 PetscFunctionReturn(0); 654 } 655 656 /*@C 657 658 PCPatchSetComputeOperator - Set the callback used to compute patch matrices 659 660 Logically collective on PC 661 662 Input Parameters: 663 + pc - The PC 664 . func - The callback 665 - ctx - The user context 666 667 Calling sequence of func: 668 $ func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx) 669 670 + pc - The PC 671 . point - The point 672 . x - The input solution (not used in linear problems) 673 . mat - The patch matrix 674 . cellIS - An array of the cell numbers 675 . n - The size of dofsArray 676 . dofsArray - The dofmap for the dofs to be solved for 677 . dofsArrayWithAll - The dofmap for all dofs on the patch 678 - ctx - The user context 679 680 Level: advanced 681 682 Notes: 683 The matrix entries have been set to zero before the call. 684 685 .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo() 686 @*/ 687 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 688 { 689 PC_PATCH *patch = (PC_PATCH *) pc->data; 690 691 PetscFunctionBegin; 692 patch->usercomputeop = func; 693 patch->usercomputeopctx = ctx; 694 PetscFunctionReturn(0); 695 } 696 697 /*@C 698 699 PCPatchSetComputeOperatorInteriorFacets - Set the callback used to compute facet integrals for patch matrices 700 701 Logically collective on PC 702 703 Input Parameters: 704 + pc - The PC 705 . func - The callback 706 - ctx - The user context 707 708 Calling sequence of func: 709 $ func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx) 710 711 + pc - The PC 712 . point - The point 713 . x - The input solution (not used in linear problems) 714 . mat - The patch matrix 715 . facetIS - An array of the facet numbers 716 . n - The size of dofsArray 717 . dofsArray - The dofmap for the dofs to be solved for 718 . dofsArrayWithAll - The dofmap for all dofs on the patch 719 - ctx - The user context 720 721 Level: advanced 722 723 Notes: 724 The matrix entries have been set to zero before the call. 725 726 .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo() 727 @*/ 728 PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 729 { 730 PC_PATCH *patch = (PC_PATCH *) pc->data; 731 732 PetscFunctionBegin; 733 patch->usercomputeopintfacet = func; 734 patch->usercomputeopintfacetctx = ctx; 735 PetscFunctionReturn(0); 736 } 737 738 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 739 on exit, cht contains all the topological entities we need to compute their residuals. 740 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 741 here we assume a standard FE sparsity pattern.*/ 742 /* TODO: Use DMPlexGetAdjacency() */ 743 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht) 744 { 745 DM dm; 746 PC_PATCH *patch = (PC_PATCH *) pc->data; 747 PetscHashIter hi; 748 PetscInt point; 749 PetscInt *star = NULL, *closure = NULL; 750 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 751 PetscInt *fStar = NULL, *fClosure = NULL; 752 PetscInt fBegin, fEnd, fsi, fci, fStarSize, fClosureSize; 753 PetscErrorCode ierr; 754 755 PetscFunctionBegin; 756 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 757 ierr = DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd);CHKERRQ(ierr); 758 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 759 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 760 ierr = PetscHSetIClear(cht);CHKERRQ(ierr); 761 PetscHashIterBegin(ht, hi); 762 while (!PetscHashIterAtEnd(ht, hi)) { 763 764 PetscHashIterGetKey(ht, hi, point); 765 PetscHashIterNext(ht, hi); 766 767 /* Loop over all the cells that this point connects to */ 768 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 769 for (si = 0; si < starSize*2; si += 2) { 770 const PetscInt ownedpoint = star[si]; 771 /* TODO Check for point in cht before running through closure again */ 772 /* now loop over all entities in the closure of that cell */ 773 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 774 for (ci = 0; ci < closureSize*2; ci += 2) { 775 const PetscInt seenpoint = closure[ci]; 776 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 777 ierr = PetscHSetIAdd(cht, seenpoint);CHKERRQ(ierr); 778 /* Facet integrals couple dofs across facets, so in that case for each of 779 * the facets we need to add all dofs on the other side of the facet to 780 * the seen dofs. */ 781 if(patch->usercomputeopintfacet){ 782 if(fBegin <= seenpoint && seenpoint < fEnd){ 783 ierr = DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar);CHKERRQ(ierr); 784 for (fsi = 0; fsi < fStarSize*2; fsi += 2) { 785 ierr = DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure);CHKERRQ(ierr); 786 for (fci = 0; fci < fClosureSize*2; fci += 2) { 787 ierr = PetscHSetIAdd(cht, fClosure[fci]);CHKERRQ(ierr); 788 } 789 ierr = DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure);CHKERRQ(ierr); 790 } 791 ierr = DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar);CHKERRQ(ierr); 792 } 793 } 794 } 795 ierr = DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 796 } 797 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 798 } 799 PetscFunctionReturn(0); 800 } 801 802 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 803 { 804 PetscErrorCode ierr; 805 806 PetscFunctionBegin; 807 if (combined) { 808 if (f < 0) { 809 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 810 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 811 } else { 812 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 813 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 814 } 815 } else { 816 if (f < 0) { 817 PC_PATCH *patch = (PC_PATCH *) pc->data; 818 PetscInt fdof, g; 819 820 if (dof) { 821 *dof = 0; 822 for (g = 0; g < patch->nsubspaces; ++g) { 823 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 824 *dof += fdof; 825 } 826 } 827 if (off) { 828 *off = 0; 829 for (g = 0; g < patch->nsubspaces; ++g) { 830 ierr = PetscSectionGetOffset(dofSection[g], p, &fdof);CHKERRQ(ierr); 831 *off += fdof; 832 } 833 } 834 } else { 835 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 836 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 837 } 838 } 839 PetscFunctionReturn(0); 840 } 841 842 /* Given a hash table with a set of topological entities (pts), compute the degrees of 843 freedom in global concatenated numbering on those entities. 844 For Vanka smoothing, this needs to do something special: ignore dofs of the 845 constraint subspace on entities that aren't the base entity we're building the patch 846 around. */ 847 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude) 848 { 849 PC_PATCH *patch = (PC_PATCH *) pc->data; 850 PetscHashIter hi; 851 PetscInt ldof, loff; 852 PetscInt k, p; 853 PetscErrorCode ierr; 854 855 PetscFunctionBegin; 856 ierr = PetscHSetIClear(dofs);CHKERRQ(ierr); 857 for (k = 0; k < patch->nsubspaces; ++k) { 858 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 859 PetscInt bs = patch->bs[k]; 860 PetscInt j, l; 861 862 if (subspaces_to_exclude != NULL) { 863 PetscBool should_exclude_k = PETSC_FALSE; 864 PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k); 865 if (should_exclude_k) { 866 /* only get this subspace dofs at the base entity, not any others */ 867 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 868 if (0 == ldof) continue; 869 for (j = loff; j < ldof + loff; ++j) { 870 for (l = 0; l < bs; ++l) { 871 PetscInt dof = bs*j + l + subspaceOffset; 872 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 873 } 874 } 875 continue; /* skip the other dofs of this subspace */ 876 } 877 } 878 879 PetscHashIterBegin(pts, hi); 880 while (!PetscHashIterAtEnd(pts, hi)) { 881 PetscHashIterGetKey(pts, hi, p); 882 PetscHashIterNext(pts, hi); 883 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 884 if (0 == ldof) continue; 885 for (j = loff; j < ldof + loff; ++j) { 886 for (l = 0; l < bs; ++l) { 887 PetscInt dof = bs*j + l + subspaceOffset; 888 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 889 } 890 } 891 } 892 } 893 PetscFunctionReturn(0); 894 } 895 896 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 897 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C) 898 { 899 PetscHashIter hi; 900 PetscInt key; 901 PetscBool flg; 902 PetscErrorCode ierr; 903 904 PetscFunctionBegin; 905 ierr = PetscHSetIClear(C);CHKERRQ(ierr); 906 PetscHashIterBegin(B, hi); 907 while (!PetscHashIterAtEnd(B, hi)) { 908 PetscHashIterGetKey(B, hi, key); 909 PetscHashIterNext(B, hi); 910 ierr = PetscHSetIHas(A, key, &flg);CHKERRQ(ierr); 911 if (!flg) {ierr = PetscHSetIAdd(C, key);CHKERRQ(ierr);} 912 } 913 PetscFunctionReturn(0); 914 } 915 916 /* 917 * PCPatchCreateCellPatches - create patches. 918 * 919 * Input Parameters: 920 * + dm - The DMPlex object defining the mesh 921 * 922 * Output Parameters: 923 * + cellCounts - Section with counts of cells around each vertex 924 * . cells - IS of the cell point indices of cells in each patch 925 * . pointCounts - Section with counts of cells around each vertex 926 * - point - IS of the cell point indices of cells in each patch 927 */ 928 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 929 { 930 PC_PATCH *patch = (PC_PATCH *) pc->data; 931 DMLabel ghost = NULL; 932 DM dm, plex; 933 PetscHSetI ht, cht; 934 PetscSection cellCounts, pointCounts, intFacetCounts, extFacetCounts; 935 PetscInt *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell; 936 PetscInt numCells, numPoints, numIntFacets, numExtFacets; 937 const PetscInt *leaves; 938 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 939 PetscBool isFiredrake; 940 PetscErrorCode ierr; 941 942 PetscFunctionBegin; 943 /* Used to keep track of the cells in the patch. */ 944 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 945 ierr = PetscHSetICreate(&cht);CHKERRQ(ierr); 946 947 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 948 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 949 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 950 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 951 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 952 953 if (patch->user_patches) { 954 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 955 vStart = 0; vEnd = patch->npatch; 956 } else if (patch->ctype == PC_PATCH_PARDECOMP) { 957 vStart = 0; vEnd = 1; 958 } else if (patch->codim < 0) { 959 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 960 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 961 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 962 patch->npatch = vEnd - vStart; 963 964 /* These labels mark the owned points. We only create patches around points that this process owns. */ 965 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 966 if (isFiredrake) { 967 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 968 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 969 } else { 970 PetscSF sf; 971 972 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 973 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 974 nleaves = PetscMax(nleaves, 0); 975 } 976 977 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 978 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 979 cellCounts = patch->cellCounts; 980 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 981 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 982 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 983 pointCounts = patch->pointCounts; 984 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 985 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts);CHKERRQ(ierr); 986 ierr = PetscObjectSetName((PetscObject) patch->extFacetCounts, "Patch Exterior Facet Layout");CHKERRQ(ierr); 987 extFacetCounts = patch->extFacetCounts; 988 ierr = PetscSectionSetChart(extFacetCounts, vStart, vEnd);CHKERRQ(ierr); 989 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts);CHKERRQ(ierr); 990 ierr = PetscObjectSetName((PetscObject) patch->intFacetCounts, "Patch Interior Facet Layout");CHKERRQ(ierr); 991 intFacetCounts = patch->intFacetCounts; 992 ierr = PetscSectionSetChart(intFacetCounts, vStart, vEnd);CHKERRQ(ierr); 993 /* Count cells and points in the patch surrounding each entity */ 994 PetscInt fStart, fEnd; 995 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 996 for (v = vStart; v < vEnd; ++v) { 997 PetscHashIter hi; 998 PetscInt chtSize, loc = -1; 999 PetscBool flg; 1000 1001 if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) { 1002 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 1003 else {ierr = PetscFindInt(v, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 1004 /* Not an owned entity, don't make a cell patch. */ 1005 if (flg) continue; 1006 } 1007 1008 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 1009 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 1010 ierr = PetscHSetIGetSize(cht, &chtSize);CHKERRQ(ierr); 1011 /* empty patch, continue */ 1012 if (chtSize == 0) continue; 1013 1014 /* safe because size(cht) > 0 from above */ 1015 PetscHashIterBegin(cht, hi); 1016 while (!PetscHashIterAtEnd(cht, hi)) { 1017 PetscInt point, pdof; 1018 1019 PetscHashIterGetKey(cht, hi, point); 1020 if (fStart <= point && point < fEnd) { 1021 const PetscInt *support; 1022 PetscInt supportSize, p; 1023 PetscBool interior = PETSC_TRUE; 1024 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1025 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1026 if (supportSize == 1) { 1027 interior = PETSC_FALSE; 1028 } else { 1029 for (p = 0; p < supportSize; p++) { 1030 PetscBool found; 1031 /* FIXME: can I do this while iterating over cht? */ 1032 PetscHSetIHas(cht, support[p], &found); 1033 if (!found) { 1034 interior = PETSC_FALSE; 1035 break; 1036 } 1037 } 1038 } 1039 if (interior) { 1040 ierr = PetscSectionAddDof(intFacetCounts, v, 1);CHKERRQ(ierr); 1041 } else { 1042 ierr = PetscSectionAddDof(extFacetCounts, v, 1);CHKERRQ(ierr); 1043 } 1044 } 1045 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 1046 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 1047 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 1048 PetscHashIterNext(cht, hi); 1049 } 1050 } 1051 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 1052 1053 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 1054 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 1055 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 1056 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 1057 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 1058 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 1059 1060 ierr = PetscSectionSetUp(intFacetCounts);CHKERRQ(ierr); 1061 ierr = PetscSectionSetUp(extFacetCounts);CHKERRQ(ierr); 1062 ierr = PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);CHKERRQ(ierr); 1063 ierr = PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);CHKERRQ(ierr); 1064 ierr = PetscMalloc1(numIntFacets, &intFacetsArray);CHKERRQ(ierr); 1065 ierr = PetscMalloc1(numIntFacets*2, &intFacetsToPatchCell);CHKERRQ(ierr); 1066 ierr = PetscMalloc1(numExtFacets, &extFacetsArray);CHKERRQ(ierr); 1067 1068 1069 /* Now that we know how much space we need, run through again and actually remember the cells. */ 1070 for (v = vStart; v < vEnd; v++ ) { 1071 PetscHashIter hi; 1072 PetscInt dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0; 1073 1074 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 1075 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 1076 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 1077 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 1078 ierr = PetscSectionGetDof(intFacetCounts, v, &ifdof);CHKERRQ(ierr); 1079 ierr = PetscSectionGetOffset(intFacetCounts, v, &ifoff);CHKERRQ(ierr); 1080 ierr = PetscSectionGetDof(extFacetCounts, v, &efdof);CHKERRQ(ierr); 1081 ierr = PetscSectionGetOffset(extFacetCounts, v, &efoff);CHKERRQ(ierr); 1082 if (dof <= 0) continue; 1083 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 1084 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 1085 PetscHashIterBegin(cht, hi); 1086 while (!PetscHashIterAtEnd(cht, hi)) { 1087 PetscInt point; 1088 1089 PetscHashIterGetKey(cht, hi, point); 1090 if (fStart <= point && point < fEnd) { 1091 const PetscInt *support; 1092 PetscInt supportSize, p; 1093 PetscBool interior = PETSC_TRUE; 1094 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1095 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1096 if (supportSize == 1) { 1097 interior = PETSC_FALSE; 1098 } else { 1099 for (p = 0; p < supportSize; p++) { 1100 PetscBool found; 1101 /* FIXME: can I do this while iterating over cht? */ 1102 PetscHSetIHas(cht, support[p], &found); 1103 if (!found) { 1104 interior = PETSC_FALSE; 1105 break; 1106 } 1107 } 1108 } 1109 if (interior) { 1110 intFacetsToPatchCell[2*(ifoff + ifn)] = support[0]; 1111 intFacetsToPatchCell[2*(ifoff + ifn) + 1] = support[1]; 1112 intFacetsArray[ifoff + ifn++] = point; 1113 } else { 1114 extFacetsArray[efoff + efn++] = point; 1115 } 1116 } 1117 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 1118 if (pdof) {pointsArray[off + n++] = point;} 1119 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 1120 PetscHashIterNext(cht, hi); 1121 } 1122 if (ifn != ifdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of interior facets in patch %D is %D, but should be %D", v, ifn, ifdof); 1123 if (efn != efdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of exterior facets in patch %D is %D, but should be %D", v, efn, efdof); 1124 if (cn != cdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %D is %D, but should be %D", v, cn, cdof); 1125 if (n != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %D is %D, but should be %D", v, n, dof); 1126 1127 for (ifn = 0; ifn < ifdof; ifn++) { 1128 PetscInt cell0 = intFacetsToPatchCell[2*(ifoff + ifn)]; 1129 PetscInt cell1 = intFacetsToPatchCell[2*(ifoff + ifn) + 1]; 1130 PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE; 1131 for (n = 0; n < cdof; n++) { 1132 if (!found0 && cell0 == cellsArray[coff + n]) { 1133 intFacetsToPatchCell[2*(ifoff + ifn)] = n; 1134 found0 = PETSC_TRUE; 1135 } 1136 if (!found1 && cell1 == cellsArray[coff + n]) { 1137 intFacetsToPatchCell[2*(ifoff + ifn) + 1] = n; 1138 found1 = PETSC_TRUE; 1139 } 1140 if (found0 && found1) break; 1141 } 1142 if (!(found0 && found1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support"); 1143 } 1144 } 1145 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 1146 ierr = PetscHSetIDestroy(&cht);CHKERRQ(ierr); 1147 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1148 1149 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 1150 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 1151 if (patch->viewCells) { 1152 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 1153 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 1154 } 1155 ierr = ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets);CHKERRQ(ierr); 1156 ierr = PetscObjectSetName((PetscObject) patch->intFacets, "Patch Interior Facets");CHKERRQ(ierr); 1157 ierr = ISCreateGeneral(PETSC_COMM_SELF, 2*numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell);CHKERRQ(ierr); 1158 ierr = PetscObjectSetName((PetscObject) patch->intFacetsToPatchCell, "Patch Interior Facets local support");CHKERRQ(ierr); 1159 if (patch->viewIntFacets) { 1160 ierr = ObjectView((PetscObject) patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets);CHKERRQ(ierr); 1161 ierr = ObjectView((PetscObject) patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets);CHKERRQ(ierr); 1162 ierr = ObjectView((PetscObject) patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets);CHKERRQ(ierr); 1163 } 1164 ierr = ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets);CHKERRQ(ierr); 1165 ierr = PetscObjectSetName((PetscObject) patch->extFacets, "Patch Exterior Facets");CHKERRQ(ierr); 1166 if (patch->viewExtFacets) { 1167 ierr = ObjectView((PetscObject) patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets);CHKERRQ(ierr); 1168 ierr = ObjectView((PetscObject) patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets);CHKERRQ(ierr); 1169 } 1170 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 1171 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 1172 if (patch->viewPoints) { 1173 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 1174 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 1175 } 1176 PetscFunctionReturn(0); 1177 } 1178 1179 /* 1180 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 1181 * 1182 * Input Parameters: 1183 * + dm - The DMPlex object defining the mesh 1184 * . cellCounts - Section with counts of cells around each vertex 1185 * . cells - IS of the cell point indices of cells in each patch 1186 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 1187 * . nodesPerCell - number of nodes per cell. 1188 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 1189 * 1190 * Output Parameters: 1191 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 1192 * . gtolCounts - Section with counts of dofs per cell patch 1193 * - gtol - IS mapping from global dofs to local dofs for each patch. 1194 */ 1195 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 1196 { 1197 PC_PATCH *patch = (PC_PATCH *) pc->data; 1198 PetscSection cellCounts = patch->cellCounts; 1199 PetscSection pointCounts = patch->pointCounts; 1200 PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL; 1201 IS cells = patch->cells; 1202 IS points = patch->points; 1203 PetscSection cellNumbering = patch->cellNumbering; 1204 PetscInt Nf = patch->nsubspaces; 1205 PetscInt numCells, numPoints; 1206 PetscInt numDofs; 1207 PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll; 1208 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 1209 PetscInt vStart, vEnd, v; 1210 const PetscInt *cellsArray, *pointsArray; 1211 PetscInt *newCellsArray = NULL; 1212 PetscInt *dofsArray = NULL; 1213 PetscInt *dofsArrayWithArtificial = NULL; 1214 PetscInt *dofsArrayWithAll = NULL; 1215 PetscInt *offsArray = NULL; 1216 PetscInt *offsArrayWithArtificial = NULL; 1217 PetscInt *offsArrayWithAll = NULL; 1218 PetscInt *asmArray = NULL; 1219 PetscInt *asmArrayWithArtificial = NULL; 1220 PetscInt *asmArrayWithAll = NULL; 1221 PetscInt *globalDofsArray = NULL; 1222 PetscInt *globalDofsArrayWithArtificial = NULL; 1223 PetscInt *globalDofsArrayWithAll = NULL; 1224 PetscInt globalIndex = 0; 1225 PetscInt key = 0; 1226 PetscInt asmKey = 0; 1227 DM dm = NULL; 1228 const PetscInt *bcNodes = NULL; 1229 PetscHMapI ht; 1230 PetscHMapI htWithArtificial; 1231 PetscHMapI htWithAll; 1232 PetscHSetI globalBcs; 1233 PetscInt numBcs; 1234 PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 1235 PetscInt pStart, pEnd, p, i; 1236 char option[PETSC_MAX_PATH_LEN]; 1237 PetscBool isNonlinear; 1238 PetscErrorCode ierr; 1239 1240 PetscFunctionBegin; 1241 1242 ierr = PCGetDM(pc, &dm); CHKERRQ(ierr); 1243 /* dofcounts section is cellcounts section * dofPerCell */ 1244 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 1245 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 1246 numDofs = numCells * totalDofsPerCell; 1247 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 1248 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 1249 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 1250 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 1251 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 1252 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 1253 gtolCounts = patch->gtolCounts; 1254 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 1255 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 1256 1257 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) 1258 { 1259 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);CHKERRQ(ierr); 1260 ierr = PetscMalloc1(numDofs, &asmArrayWithArtificial);CHKERRQ(ierr); 1261 ierr = PetscMalloc1(numDofs, &dofsArrayWithArtificial);CHKERRQ(ierr); 1262 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);CHKERRQ(ierr); 1263 gtolCountsWithArtificial = patch->gtolCountsWithArtificial; 1264 ierr = PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);CHKERRQ(ierr); 1265 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");CHKERRQ(ierr); 1266 } 1267 1268 isNonlinear = patch->isNonlinear; 1269 if(isNonlinear) 1270 { 1271 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithAll);CHKERRQ(ierr); 1272 ierr = PetscMalloc1(numDofs, &asmArrayWithAll);CHKERRQ(ierr); 1273 ierr = PetscMalloc1(numDofs, &dofsArrayWithAll);CHKERRQ(ierr); 1274 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);CHKERRQ(ierr); 1275 gtolCountsWithAll = patch->gtolCountsWithAll; 1276 ierr = PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);CHKERRQ(ierr); 1277 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");CHKERRQ(ierr); 1278 } 1279 1280 /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet 1281 conditions */ 1282 ierr = PetscHSetICreate(&globalBcs);CHKERRQ(ierr); 1283 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 1284 ierr = ISGetSize(patch->ghostBcNodes, &numBcs); CHKERRQ(ierr); 1285 for (i = 0; i < numBcs; ++i) { 1286 ierr = PetscHSetIAdd(globalBcs, bcNodes[i]);CHKERRQ(ierr); /* these are already in concatenated numbering */ 1287 } 1288 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 1289 ierr = ISDestroy(&patch->ghostBcNodes); CHKERRQ(ierr); /* memory optimisation */ 1290 1291 /* Hash tables for artificial BC construction */ 1292 ierr = PetscHSetICreate(&ownedpts);CHKERRQ(ierr); 1293 ierr = PetscHSetICreate(&seenpts);CHKERRQ(ierr); 1294 ierr = PetscHSetICreate(&owneddofs);CHKERRQ(ierr); 1295 ierr = PetscHSetICreate(&seendofs);CHKERRQ(ierr); 1296 ierr = PetscHSetICreate(&artificialbcs);CHKERRQ(ierr); 1297 1298 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 1299 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 1300 ierr = PetscHMapICreate(&ht);CHKERRQ(ierr); 1301 ierr = PetscHMapICreate(&htWithArtificial);CHKERRQ(ierr); 1302 ierr = PetscHMapICreate(&htWithAll);CHKERRQ(ierr); 1303 for (v = vStart; v < vEnd; ++v) { 1304 PetscInt localIndex = 0; 1305 PetscInt localIndexWithArtificial = 0; 1306 PetscInt localIndexWithAll = 0; 1307 PetscInt dof, off, i, j, k, l; 1308 1309 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1310 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 1311 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 1312 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1313 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1314 if (dof <= 0) continue; 1315 1316 /* Calculate the global numbers of the artificial BC dofs here first */ 1317 ierr = patch->patchconstructop((void*)patch, dm, v, ownedpts); CHKERRQ(ierr); 1318 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts); CHKERRQ(ierr); 1319 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); CHKERRQ(ierr); 1320 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); CHKERRQ(ierr); 1321 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); CHKERRQ(ierr); 1322 if (patch->viewPatches) { 1323 PetscHSetI globalbcdofs; 1324 PetscHashIter hi; 1325 MPI_Comm comm = PetscObjectComm((PetscObject)pc); 1326 1327 ierr = PetscHSetICreate(&globalbcdofs);CHKERRQ(ierr); 1328 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); CHKERRQ(ierr); 1329 PetscHashIterBegin(owneddofs, hi); 1330 while (!PetscHashIterAtEnd(owneddofs, hi)) { 1331 PetscInt globalDof; 1332 1333 PetscHashIterGetKey(owneddofs, hi, globalDof); 1334 PetscHashIterNext(owneddofs, hi); 1335 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1336 } 1337 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1338 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); CHKERRQ(ierr); 1339 PetscHashIterBegin(seendofs, hi); 1340 while (!PetscHashIterAtEnd(seendofs, hi)) { 1341 PetscInt globalDof; 1342 PetscBool flg; 1343 1344 PetscHashIterGetKey(seendofs, hi, globalDof); 1345 PetscHashIterNext(seendofs, hi); 1346 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1347 1348 ierr = PetscHSetIHas(globalBcs, globalDof, &flg);CHKERRQ(ierr); 1349 if (flg) {ierr = PetscHSetIAdd(globalbcdofs, globalDof);CHKERRQ(ierr);} 1350 } 1351 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1352 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); CHKERRQ(ierr); 1353 ierr = PetscHSetIGetSize(globalbcdofs, &numBcs);CHKERRQ(ierr); 1354 if (numBcs > 0) { 1355 PetscHashIterBegin(globalbcdofs, hi); 1356 while (!PetscHashIterAtEnd(globalbcdofs, hi)) { 1357 PetscInt globalDof; 1358 PetscHashIterGetKey(globalbcdofs, hi, globalDof); 1359 PetscHashIterNext(globalbcdofs, hi); 1360 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1361 } 1362 } 1363 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1364 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 1365 ierr = PetscHSetIGetSize(artificialbcs, &numBcs);CHKERRQ(ierr); 1366 if (numBcs > 0) { 1367 PetscHashIterBegin(artificialbcs, hi); 1368 while (!PetscHashIterAtEnd(artificialbcs, hi)) { 1369 PetscInt globalDof; 1370 PetscHashIterGetKey(artificialbcs, hi, globalDof); 1371 PetscHashIterNext(artificialbcs, hi); 1372 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1373 } 1374 } 1375 ierr = PetscSynchronizedPrintf(comm, "\n\n"); CHKERRQ(ierr); 1376 ierr = PetscHSetIDestroy(&globalbcdofs);CHKERRQ(ierr); 1377 } 1378 for (k = 0; k < patch->nsubspaces; ++k) { 1379 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1380 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1381 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1382 PetscInt bs = patch->bs[k]; 1383 1384 for (i = off; i < off + dof; ++i) { 1385 /* Walk over the cells in this patch. */ 1386 const PetscInt c = cellsArray[i]; 1387 PetscInt cell = c; 1388 1389 /* TODO Change this to an IS */ 1390 if (cellNumbering) { 1391 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 1392 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 1393 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 1394 } 1395 newCellsArray[i] = cell; 1396 for (j = 0; j < nodesPerCell; ++j) { 1397 /* For each global dof, map it into contiguous local storage. */ 1398 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 1399 /* finally, loop over block size */ 1400 for (l = 0; l < bs; ++l) { 1401 PetscInt localDof; 1402 PetscBool isGlobalBcDof, isArtificialBcDof; 1403 1404 /* first, check if this is either a globally enforced or locally enforced BC dof */ 1405 ierr = PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);CHKERRQ(ierr); 1406 ierr = PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);CHKERRQ(ierr); 1407 1408 /* if it's either, don't ever give it a local dof number */ 1409 if (isGlobalBcDof || isArtificialBcDof) { 1410 dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */ 1411 } else { 1412 ierr = PetscHMapIGet(ht, globalDof + l, &localDof);CHKERRQ(ierr); 1413 if (localDof == -1) { 1414 localDof = localIndex++; 1415 ierr = PetscHMapISet(ht, globalDof + l, localDof);CHKERRQ(ierr); 1416 } 1417 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1418 /* And store. */ 1419 dofsArray[globalIndex] = localDof; 1420 } 1421 1422 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1423 if (isGlobalBcDof) { 1424 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */ 1425 } else { 1426 ierr = PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);CHKERRQ(ierr); 1427 if (localDof == -1) { 1428 localDof = localIndexWithArtificial++; 1429 ierr = PetscHMapISet(htWithArtificial, globalDof + l, localDof);CHKERRQ(ierr); 1430 } 1431 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1432 /* And store.*/ 1433 dofsArrayWithArtificial[globalIndex] = localDof; 1434 } 1435 } 1436 1437 if(isNonlinear) { 1438 /* Build the dofmap for the function space with _all_ dofs, 1439 including those in any kind of boundary condition */ 1440 ierr = PetscHMapIGet(htWithAll, globalDof + l, &localDof);CHKERRQ(ierr); 1441 if (localDof == -1) { 1442 localDof = localIndexWithAll++; 1443 ierr = PetscHMapISet(htWithAll, globalDof + l, localDof);CHKERRQ(ierr); 1444 } 1445 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1446 /* And store.*/ 1447 dofsArrayWithAll[globalIndex] = localDof; 1448 } 1449 globalIndex++; 1450 } 1451 } 1452 } 1453 } 1454 /*How many local dofs in this patch? */ 1455 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1456 ierr = PetscHMapIGetSize(htWithArtificial, &dof);CHKERRQ(ierr); 1457 ierr = PetscSectionSetDof(gtolCountsWithArtificial, v, dof);CHKERRQ(ierr); 1458 } 1459 if (isNonlinear) { 1460 ierr = PetscHMapIGetSize(htWithAll, &dof);CHKERRQ(ierr); 1461 ierr = PetscSectionSetDof(gtolCountsWithAll, v, dof);CHKERRQ(ierr); 1462 } 1463 ierr = PetscHMapIGetSize(ht, &dof);CHKERRQ(ierr); 1464 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 1465 } 1466 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 1467 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 1468 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 1469 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 1470 1471 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1472 ierr = PetscSectionSetUp(gtolCountsWithArtificial);CHKERRQ(ierr); 1473 ierr = PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);CHKERRQ(ierr); 1474 ierr = PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);CHKERRQ(ierr); 1475 } 1476 if (isNonlinear) { 1477 ierr = PetscSectionSetUp(gtolCountsWithAll);CHKERRQ(ierr); 1478 ierr = PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);CHKERRQ(ierr); 1479 ierr = PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);CHKERRQ(ierr); 1480 } 1481 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 1482 for (v = vStart; v < vEnd; ++v) { 1483 PetscHashIter hi; 1484 PetscInt dof, off, Np, ooff, i, j, k, l; 1485 1486 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1487 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 1488 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 1489 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1490 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1491 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 1492 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 1493 if (dof <= 0) continue; 1494 1495 for (k = 0; k < patch->nsubspaces; ++k) { 1496 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1497 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1498 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1499 PetscInt bs = patch->bs[k]; 1500 PetscInt goff; 1501 1502 for (i = off; i < off + dof; ++i) { 1503 /* Reconstruct mapping of global-to-local on this patch. */ 1504 const PetscInt c = cellsArray[i]; 1505 PetscInt cell = c; 1506 1507 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1508 for (j = 0; j < nodesPerCell; ++j) { 1509 for (l = 0; l < bs; ++l) { 1510 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1511 const PetscInt localDof = dofsArray[key]; 1512 if (localDof >= 0) {ierr = PetscHMapISet(ht, globalDof, localDof);CHKERRQ(ierr);} 1513 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1514 const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key]; 1515 if (localDofWithArtificial >= 0) { 1516 ierr = PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);CHKERRQ(ierr); 1517 } 1518 } 1519 if (isNonlinear) { 1520 const PetscInt localDofWithAll = dofsArrayWithAll[key]; 1521 if (localDofWithAll >= 0) { 1522 ierr = PetscHMapISet(htWithAll, globalDof, localDofWithAll);CHKERRQ(ierr); 1523 } 1524 } 1525 key++; 1526 } 1527 } 1528 } 1529 1530 /* Shove it in the output data structure. */ 1531 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 1532 PetscHashIterBegin(ht, hi); 1533 while (!PetscHashIterAtEnd(ht, hi)) { 1534 PetscInt globalDof, localDof; 1535 1536 PetscHashIterGetKey(ht, hi, globalDof); 1537 PetscHashIterGetVal(ht, hi, localDof); 1538 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 1539 PetscHashIterNext(ht, hi); 1540 } 1541 1542 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1543 ierr = PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);CHKERRQ(ierr); 1544 PetscHashIterBegin(htWithArtificial, hi); 1545 while (!PetscHashIterAtEnd(htWithArtificial, hi)) { 1546 PetscInt globalDof, localDof; 1547 PetscHashIterGetKey(htWithArtificial, hi, globalDof); 1548 PetscHashIterGetVal(htWithArtificial, hi, localDof); 1549 if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof; 1550 PetscHashIterNext(htWithArtificial, hi); 1551 } 1552 } 1553 if (isNonlinear) { 1554 ierr = PetscSectionGetOffset(gtolCountsWithAll, v, &goff);CHKERRQ(ierr); 1555 PetscHashIterBegin(htWithAll, hi); 1556 while (!PetscHashIterAtEnd(htWithAll, hi)) { 1557 PetscInt globalDof, localDof; 1558 PetscHashIterGetKey(htWithAll, hi, globalDof); 1559 PetscHashIterGetVal(htWithAll, hi, localDof); 1560 if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof; 1561 PetscHashIterNext(htWithAll, hi); 1562 } 1563 } 1564 1565 for (p = 0; p < Np; ++p) { 1566 const PetscInt point = pointsArray[ooff + p]; 1567 PetscInt globalDof, localDof; 1568 1569 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 1570 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1571 offsArray[(ooff + p)*Nf + k] = localDof; 1572 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1573 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1574 offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof; 1575 } 1576 if (isNonlinear) { 1577 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1578 offsArrayWithAll[(ooff + p)*Nf + k] = localDof; 1579 } 1580 } 1581 } 1582 1583 ierr = PetscHSetIDestroy(&globalBcs);CHKERRQ(ierr); 1584 ierr = PetscHSetIDestroy(&ownedpts);CHKERRQ(ierr); 1585 ierr = PetscHSetIDestroy(&seenpts);CHKERRQ(ierr); 1586 ierr = PetscHSetIDestroy(&owneddofs);CHKERRQ(ierr); 1587 ierr = PetscHSetIDestroy(&seendofs);CHKERRQ(ierr); 1588 ierr = PetscHSetIDestroy(&artificialbcs);CHKERRQ(ierr); 1589 1590 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 1591 We need to create the dof table laid out cellwise first, then by subspace, 1592 as the assembler assembles cell-wise and we need to stuff the different 1593 contributions of the different function spaces to the right places. So we loop 1594 over cells, then over subspaces. */ 1595 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 1596 for (i = off; i < off + dof; ++i) { 1597 const PetscInt c = cellsArray[i]; 1598 PetscInt cell = c; 1599 1600 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1601 for (k = 0; k < patch->nsubspaces; ++k) { 1602 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1603 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1604 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1605 PetscInt bs = patch->bs[k]; 1606 1607 for (j = 0; j < nodesPerCell; ++j) { 1608 for (l = 0; l < bs; ++l) { 1609 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1610 PetscInt localDof; 1611 1612 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1613 /* If it's not in the hash table, i.e. is a BC dof, 1614 then the PetscHSetIMap above gives -1, which matches 1615 exactly the convention for PETSc's matrix assembly to 1616 ignore the dof. So we don't need to do anything here */ 1617 asmArray[asmKey] = localDof; 1618 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1619 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1620 asmArrayWithArtificial[asmKey] = localDof; 1621 } 1622 if (isNonlinear) { 1623 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1624 asmArrayWithAll[asmKey] = localDof; 1625 } 1626 asmKey++; 1627 } 1628 } 1629 } 1630 } 1631 } 1632 } 1633 if (1 == patch->nsubspaces) { 1634 ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1635 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1636 ierr = PetscMemcpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1637 } 1638 if (isNonlinear) { 1639 ierr = PetscMemcpy(asmArrayWithAll, dofsArrayWithAll, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1640 } 1641 } 1642 1643 ierr = PetscHMapIDestroy(&ht);CHKERRQ(ierr); 1644 ierr = PetscHMapIDestroy(&htWithArtificial);CHKERRQ(ierr); 1645 ierr = PetscHMapIDestroy(&htWithAll);CHKERRQ(ierr); 1646 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1647 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1648 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1649 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1650 ierr = PetscFree(dofsArrayWithArtificial);CHKERRQ(ierr); 1651 } 1652 if (isNonlinear) { 1653 ierr = PetscFree(dofsArrayWithAll);CHKERRQ(ierr); 1654 } 1655 /* Create placeholder section for map from points to patch dofs */ 1656 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1657 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1658 if (patch->combined) { 1659 PetscInt numFields; 1660 ierr = PetscSectionGetNumFields(patch->dofSection[0], &numFields);CHKERRQ(ierr); 1661 if (numFields != patch->nsubspaces) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %D and number of subspaces %D", numFields, patch->nsubspaces); 1662 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1663 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1664 for (p = pStart; p < pEnd; ++p) { 1665 PetscInt dof, fdof, f; 1666 1667 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1668 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1669 for (f = 0; f < patch->nsubspaces; ++f) { 1670 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1671 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1672 } 1673 } 1674 } else { 1675 PetscInt pStartf, pEndf, f; 1676 pStart = PETSC_MAX_INT; 1677 pEnd = PETSC_MIN_INT; 1678 for (f = 0; f < patch->nsubspaces; ++f) { 1679 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1680 pStart = PetscMin(pStart, pStartf); 1681 pEnd = PetscMax(pEnd, pEndf); 1682 } 1683 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1684 for (f = 0; f < patch->nsubspaces; ++f) { 1685 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1686 for (p = pStartf; p < pEndf; ++p) { 1687 PetscInt fdof; 1688 ierr = PetscSectionGetDof(patch->dofSection[f], p, &fdof);CHKERRQ(ierr); 1689 ierr = PetscSectionAddDof(patch->patchSection, p, fdof);CHKERRQ(ierr); 1690 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1691 } 1692 } 1693 } 1694 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1695 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1696 /* Replace cell indices with firedrake-numbered ones. */ 1697 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1698 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1699 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1700 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);CHKERRQ(ierr); 1701 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);CHKERRQ(ierr); 1702 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, option);CHKERRQ(ierr); 1703 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1704 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1705 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1706 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);CHKERRQ(ierr); 1707 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);CHKERRQ(ierr); 1708 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);CHKERRQ(ierr); 1709 } 1710 if (isNonlinear) { 1711 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);CHKERRQ(ierr); 1712 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);CHKERRQ(ierr); 1713 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);CHKERRQ(ierr); 1714 } 1715 PetscFunctionReturn(0); 1716 } 1717 1718 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial) 1719 { 1720 PC_PATCH *patch = (PC_PATCH *) pc->data; 1721 Vec x, y; 1722 PetscBool flg; 1723 PetscInt csize, rsize; 1724 const char *prefix = NULL; 1725 PetscErrorCode ierr; 1726 1727 PetscFunctionBegin; 1728 if(withArtificial) { 1729 /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */ 1730 x = patch->patchRHSWithArtificial[point]; 1731 y = patch->patchRHSWithArtificial[point]; 1732 } else { 1733 x = patch->patchRHS[point]; 1734 y = patch->patchUpdate[point]; 1735 } 1736 1737 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1738 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1739 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1740 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1741 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1742 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1743 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1744 else if (!patch->sub_mat_type) {ierr = MatSetType(*mat, MATDENSE);CHKERRQ(ierr);} 1745 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1746 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1747 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1748 /* Sparse patch matrices */ 1749 if (!flg) { 1750 PetscBT bt; 1751 PetscInt *dnnz = NULL; 1752 const PetscInt *dofsArray = NULL; 1753 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1754 1755 if(withArtificial) { 1756 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1757 } else { 1758 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1759 } 1760 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1761 point += pStart; 1762 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1763 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1764 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1765 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1766 /* A PetscBT uses N^2 bits to store the sparsity pattern on a 1767 * patch. This is probably OK if the patches are not too big, 1768 * but uses too much memory. We therefore switch based on rsize. */ 1769 if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */ 1770 PetscScalar *zeroes; 1771 PetscInt rows; 1772 1773 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1774 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1775 for (c = 0; c < ncell; ++c) { 1776 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1777 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1778 const PetscInt row = idx[i]; 1779 if (row < 0) continue; 1780 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1781 const PetscInt col = idx[j]; 1782 const PetscInt key = row*rsize + col; 1783 if (col < 0) continue; 1784 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1785 } 1786 } 1787 } 1788 1789 if (patch->usercomputeopintfacet) { 1790 const PetscInt *intFacetsArray = NULL; 1791 PetscInt i, numIntFacets, intFacetOffset; 1792 const PetscInt *facetCells = NULL; 1793 1794 ierr = PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);CHKERRQ(ierr); 1795 ierr = PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);CHKERRQ(ierr); 1796 ierr = ISGetIndices(patch->intFacetsToPatchCell, &facetCells);CHKERRQ(ierr); 1797 ierr = ISGetIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 1798 for (i = 0; i < numIntFacets; i++) { 1799 const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0]; 1800 const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1]; 1801 PetscInt celli, cellj; 1802 1803 for (celli = 0; celli < patch->totalDofsPerCell; celli++) { 1804 const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli]; 1805 if (row < 0) continue; 1806 for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) { 1807 const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj]; 1808 const PetscInt key = row*rsize + col; 1809 if (col < 0) continue; 1810 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1811 } 1812 } 1813 1814 for (celli = 0; celli < patch->totalDofsPerCell; celli++) { 1815 const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli]; 1816 if (row < 0) continue; 1817 for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) { 1818 const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj]; 1819 const PetscInt key = row*rsize + col; 1820 if (col < 0) continue; 1821 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1822 } 1823 } 1824 } 1825 } 1826 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1827 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1828 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1829 1830 ierr = PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);CHKERRQ(ierr); 1831 for (c = 0; c < ncell; ++c) { 1832 const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell]; 1833 ierr = MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);CHKERRQ(ierr); 1834 } 1835 ierr = MatGetLocalSize(*mat, &rows, NULL);CHKERRQ(ierr); 1836 for (i = 0; i < rows; ++i) { 1837 ierr = MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);CHKERRQ(ierr); 1838 } 1839 1840 if (patch->usercomputeopintfacet) { 1841 const PetscInt *intFacetsArray = NULL; 1842 PetscInt i, numIntFacets, intFacetOffset; 1843 const PetscInt *facetCells = NULL; 1844 1845 ierr = PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);CHKERRQ(ierr); 1846 ierr = PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);CHKERRQ(ierr); 1847 ierr = ISGetIndices(patch->intFacetsToPatchCell, &facetCells);CHKERRQ(ierr); 1848 ierr = ISGetIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 1849 for (i = 0; i < numIntFacets; i++) { 1850 const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0]; 1851 const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1]; 1852 const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell]; 1853 const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell]; 1854 ierr = MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);CHKERRQ(ierr); 1855 ierr = MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);CHKERRQ(ierr); 1856 } 1857 } 1858 1859 ierr = MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1860 ierr = MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1861 1862 ierr = PetscFree(zeroes);CHKERRQ(ierr); 1863 1864 } else { /* rsize too big, use MATPREALLOCATOR */ 1865 Mat preallocator; 1866 PetscScalar* vals; 1867 1868 ierr = PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);CHKERRQ(ierr); 1869 ierr = MatCreate(PETSC_COMM_SELF, &preallocator);CHKERRQ(ierr); 1870 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1871 ierr = MatSetSizes(preallocator, rsize, rsize, rsize, rsize);CHKERRQ(ierr); 1872 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1873 1874 for (c = 0; c < ncell; ++c) { 1875 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1876 ierr = MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);CHKERRQ(ierr); 1877 } 1878 1879 if (patch->usercomputeopintfacet) { 1880 const PetscInt *intFacetsArray = NULL; 1881 PetscInt i, numIntFacets, intFacetOffset; 1882 const PetscInt *facetCells = NULL; 1883 1884 ierr = PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);CHKERRQ(ierr); 1885 ierr = PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);CHKERRQ(ierr); 1886 ierr = ISGetIndices(patch->intFacetsToPatchCell, &facetCells);CHKERRQ(ierr); 1887 ierr = ISGetIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 1888 for (i = 0; i < numIntFacets; i++) { 1889 const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0]; 1890 const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1]; 1891 const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell]; 1892 const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell]; 1893 ierr = MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);CHKERRQ(ierr); 1894 ierr = MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);CHKERRQ(ierr); 1895 } 1896 } 1897 1898 ierr = PetscFree(vals);CHKERRQ(ierr); 1899 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1900 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1901 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);CHKERRQ(ierr); 1902 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1903 } 1904 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1905 if(withArtificial) { 1906 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1907 } else { 1908 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1909 } 1910 } 1911 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 1916 { 1917 PC_PATCH *patch = (PC_PATCH *) pc->data; 1918 DM dm; 1919 PetscSection s; 1920 const PetscInt *parray, *oarray; 1921 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1922 PetscErrorCode ierr; 1923 1924 PetscFunctionBegin; 1925 if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator"); 1926 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1927 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1928 /* Set offset into patch */ 1929 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1930 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1931 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1932 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1933 for (f = 0; f < Nf; ++f) { 1934 for (p = 0; p < Np; ++p) { 1935 const PetscInt point = parray[poff+p]; 1936 PetscInt dof; 1937 1938 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1939 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1940 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1941 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1942 } 1943 } 1944 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1945 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1946 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1947 ierr = DMPlexComputeResidual_Patch_Internal(pc->dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);CHKERRQ(ierr); 1948 PetscFunctionReturn(0); 1949 } 1950 1951 PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point) 1952 { 1953 PC_PATCH *patch = (PC_PATCH *) pc->data; 1954 const PetscInt *dofsArray; 1955 const PetscInt *dofsArrayWithAll; 1956 const PetscInt *cellsArray; 1957 PetscInt ncell, offset, pStart, pEnd; 1958 PetscErrorCode ierr; 1959 1960 PetscFunctionBegin; 1961 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1962 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1963 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1964 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1965 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1966 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1967 1968 point += pStart; 1969 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1970 1971 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1972 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1973 if (ncell <= 0) { 1974 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1975 PetscFunctionReturn(0); 1976 } 1977 PetscStackPush("PCPatch user callback"); 1978 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1979 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1980 ierr = patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, 1981 dofsArrayWithAll + offset*patch->totalDofsPerCell, 1982 patch->usercomputefctx);CHKERRQ(ierr); 1983 PetscStackPop; 1984 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1985 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1986 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1987 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1988 if (patch->viewMatrix) { 1989 char name[PETSC_MAX_PATH_LEN]; 1990 1991 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);CHKERRQ(ierr); 1992 ierr = PetscObjectSetName((PetscObject) F, name);CHKERRQ(ierr); 1993 ierr = ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1994 } 1995 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1996 PetscFunctionReturn(0); 1997 } 1998 1999 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 2000 { 2001 PC_PATCH *patch = (PC_PATCH *) pc->data; 2002 DM dm; 2003 PetscSection s; 2004 const PetscInt *parray, *oarray; 2005 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 2006 PetscErrorCode ierr; 2007 2008 PetscFunctionBegin; 2009 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 2010 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 2011 /* Set offset into patch */ 2012 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 2013 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 2014 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 2015 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 2016 for (f = 0; f < Nf; ++f) { 2017 for (p = 0; p < Np; ++p) { 2018 const PetscInt point = parray[poff+p]; 2019 PetscInt dof; 2020 2021 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 2022 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 2023 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 2024 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 2025 } 2026 } 2027 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 2028 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 2029 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 2030 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 2031 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);CHKERRQ(ierr); 2032 PetscFunctionReturn(0); 2033 } 2034 2035 PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial) 2036 { 2037 PC_PATCH *patch = (PC_PATCH *) pc->data; 2038 const PetscInt *dofsArray; 2039 const PetscInt *dofsArrayWithAll = NULL; 2040 const PetscInt *cellsArray; 2041 PetscInt ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset; 2042 PetscBool isNonlinear; 2043 PetscErrorCode ierr; 2044 2045 PetscFunctionBegin; 2046 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 2047 isNonlinear = patch->isNonlinear; 2048 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 2049 if(withArtificial) { 2050 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 2051 } else { 2052 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 2053 } 2054 if (isNonlinear) { 2055 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 2056 } 2057 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 2058 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 2059 2060 point += pStart; 2061 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 2062 2063 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 2064 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 2065 if (ncell <= 0) { 2066 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 2067 PetscFunctionReturn(0); 2068 } 2069 if (patch->precomputeElementTensors) { 2070 PetscInt i; 2071 PetscInt ndof = patch->totalDofsPerCell; 2072 const PetscScalar *elementTensors; 2073 2074 ierr = VecGetArrayRead(patch->cellMats, &elementTensors);CHKERRQ(ierr); 2075 for (i = 0; i < ncell; i++) { 2076 const PetscInt cell = cellsArray[i + offset]; 2077 const PetscInt *idx = dofsArray + (offset + i)*ndof; 2078 const PetscScalar *v = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof; 2079 ierr = MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);CHKERRQ(ierr); 2080 } 2081 ierr = VecRestoreArrayRead(patch->cellMats, &elementTensors);CHKERRQ(ierr); 2082 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2083 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2084 } else { 2085 PetscStackPush("PCPatch user callback"); 2086 /* Cannot reuse the same IS because the geometry info is being cached in it */ 2087 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 2088 ierr = patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);CHKERRQ(ierr); 2089 } 2090 if (patch->usercomputeopintfacet) { 2091 ierr = PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);CHKERRQ(ierr); 2092 ierr = PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);CHKERRQ(ierr); 2093 if (numIntFacets > 0) { 2094 /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */ 2095 PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL; 2096 const PetscInt *intFacetsArray = NULL; 2097 PetscInt idx = 0; 2098 PetscInt i, c, d; 2099 PetscInt fStart; 2100 DM dm; 2101 IS facetIS = NULL; 2102 const PetscInt *facetCells = NULL; 2103 2104 ierr = ISGetIndices(patch->intFacetsToPatchCell, &facetCells);CHKERRQ(ierr); 2105 ierr = ISGetIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 2106 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 2107 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, NULL);CHKERRQ(ierr); 2108 /* FIXME: Pull this malloc out. */ 2109 ierr = PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);CHKERRQ(ierr); 2110 if (dofsArrayWithAll) { 2111 ierr = PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);CHKERRQ(ierr); 2112 } 2113 if (patch->precomputeElementTensors) { 2114 PetscInt nFacetDof = 2*patch->totalDofsPerCell; 2115 const PetscScalar *elementTensors; 2116 2117 ierr = VecGetArrayRead(patch->intFacetMats, &elementTensors);CHKERRQ(ierr); 2118 2119 for (i = 0; i < numIntFacets; i++) { 2120 const PetscInt facet = intFacetsArray[i + intFacetOffset]; 2121 const PetscScalar *v = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof; 2122 idx = 0; 2123 /* 2124 * 0--1 2125 * |\-| 2126 * |+\| 2127 * 2--3 2128 * [0, 2, 3, 0, 1, 3] 2129 */ 2130 for (c = 0; c < 2; c++) { 2131 const PetscInt cell = facetCells[2*(intFacetOffset + i) + c]; 2132 for (d = 0; d < patch->totalDofsPerCell; d++) { 2133 facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d]; 2134 idx++; 2135 } 2136 } 2137 ierr = MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);CHKERRQ(ierr); 2138 } 2139 ierr = VecRestoreArrayRead(patch->intFacetMats, &elementTensors);CHKERRQ(ierr); 2140 } else { 2141 /* 2142 * 0--1 2143 * |\-| 2144 * |+\| 2145 * 2--3 2146 * [0, 2, 3, 0, 1, 3] 2147 */ 2148 for (i = 0; i < numIntFacets; i++) { 2149 for (c = 0; c < 2; c++) { 2150 const PetscInt cell = facetCells[2*(intFacetOffset + i) + c]; 2151 for (d = 0; d < patch->totalDofsPerCell; d++) { 2152 facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d]; 2153 if (dofsArrayWithAll) { 2154 facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d]; 2155 } 2156 idx++; 2157 } 2158 } 2159 } 2160 ierr = ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);CHKERRQ(ierr); 2161 ierr = patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);CHKERRQ(ierr); 2162 ierr = ISDestroy(&facetIS);CHKERRQ(ierr); 2163 } 2164 ierr = ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);CHKERRQ(ierr); 2165 ierr = ISRestoreIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 2166 ierr = PetscFree(facetDofs);CHKERRQ(ierr); 2167 ierr = PetscFree(facetDofsWithAll);CHKERRQ(ierr); 2168 } 2169 } 2170 2171 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2172 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2173 2174 PetscStackPop; 2175 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 2176 if(withArtificial) { 2177 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 2178 } else { 2179 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 2180 } 2181 if (isNonlinear) { 2182 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 2183 } 2184 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 2185 if (patch->viewMatrix) { 2186 char name[PETSC_MAX_PATH_LEN]; 2187 2188 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);CHKERRQ(ierr); 2189 ierr = PetscObjectSetName((PetscObject) mat, name);CHKERRQ(ierr); 2190 ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 2191 } 2192 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], 2197 PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv) 2198 { 2199 Vec data; 2200 PetscScalar *array; 2201 PetscInt bs, nz, i, j, cell; 2202 PetscErrorCode ierr; 2203 2204 ierr = MatShellGetContext(mat, &data);CHKERRQ(ierr); 2205 ierr = VecGetBlockSize(data, &bs);CHKERRQ(ierr); 2206 ierr = VecGetSize(data, &nz);CHKERRQ(ierr); 2207 ierr = VecGetArray(data, &array);CHKERRQ(ierr); 2208 if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion"); 2209 cell = (PetscInt)(idxm[0]/bs); // use the fact that this is called once per cell 2210 for (i = 0; i < m; i++) { 2211 //const PetscInt row = idxm[i]; 2212 if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!"); 2213 for (j = 0; j < n; j++) { 2214 const PetscScalar v_ = v[i*bs + j]; 2215 /* Indexing is special to the data structure we have! */ 2216 if (addv == INSERT_VALUES) { 2217 array[cell*bs*bs + i*bs + j] = v_; 2218 } else { 2219 array[cell*bs*bs + i*bs + j] += v_; 2220 } 2221 } 2222 } 2223 ierr = VecRestoreArray(data, &array);CHKERRQ(ierr); 2224 PetscFunctionReturn(0); 2225 } 2226 2227 static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc) 2228 { 2229 PC_PATCH *patch = (PC_PATCH *)pc->data; 2230 const PetscInt *cellsArray; 2231 PetscInt ncell, offset; 2232 const PetscInt *dofMapArray; 2233 PetscInt i, j; 2234 IS dofMap; 2235 IS cellIS; 2236 const PetscInt ndof = patch->totalDofsPerCell; 2237 PetscErrorCode ierr; 2238 Mat vecMat; 2239 PetscInt cStart, cEnd; 2240 DM dm, plex; 2241 2242 2243 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 2244 2245 if (!patch->allCells) { 2246 PetscHSetI cells; 2247 PetscHashIter hi; 2248 PetscInt pStart, pEnd; 2249 PetscInt *allCells = NULL; 2250 ierr = PetscHSetICreate(&cells);CHKERRQ(ierr); 2251 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 2252 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 2253 for (i = pStart; i < pEnd; i++) { 2254 ierr = PetscSectionGetDof(patch->cellCounts, i, &ncell);CHKERRQ(ierr); 2255 ierr = PetscSectionGetOffset(patch->cellCounts, i, &offset);CHKERRQ(ierr); 2256 if (ncell <= 0) continue; 2257 for (j = 0; j < ncell; j++) { 2258 PetscHSetIAdd(cells, cellsArray[offset + j]);CHKERRQ(ierr); 2259 } 2260 } 2261 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 2262 ierr = PetscHSetIGetSize(cells, &ncell);CHKERRQ(ierr); 2263 ierr = PetscMalloc1(ncell, &allCells);CHKERRQ(ierr); 2264 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 2265 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2266 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2267 ierr = PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);CHKERRQ(ierr); 2268 i = 0; 2269 PetscHashIterBegin(cells, hi); 2270 while (!PetscHashIterAtEnd(cells, hi)) { 2271 PetscHashIterGetKey(cells, hi, allCells[i]); 2272 patch->precomputedTensorLocations[allCells[i]] = i; 2273 PetscHashIterNext(cells, hi); 2274 i++; 2275 } 2276 ierr = PetscHSetIDestroy(&cells);CHKERRQ(ierr); 2277 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);CHKERRQ(ierr); 2278 } 2279 ierr = ISGetSize(patch->allCells, &ncell);CHKERRQ(ierr); 2280 if (!patch->cellMats) { 2281 ierr = VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);CHKERRQ(ierr); 2282 ierr = VecSetBlockSize(patch->cellMats, ndof);CHKERRQ(ierr); 2283 } 2284 ierr = VecSet(patch->cellMats, 0);CHKERRQ(ierr); 2285 2286 ierr = MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof, 2287 (void*)patch->cellMats, &vecMat); 2288 ierr = MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);CHKERRQ(ierr); 2289 ierr = ISGetSize(patch->allCells, &ncell);CHKERRQ(ierr); 2290 ierr = ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);CHKERRQ(ierr); 2291 ierr = ISGetIndices(dofMap, &dofMapArray);CHKERRQ(ierr); 2292 ierr = ISGetIndices(patch->allCells, &cellsArray);CHKERRQ(ierr); 2293 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);CHKERRQ(ierr); 2294 PetscStackPush("PCPatch user callback"); 2295 /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */ 2296 ierr = patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);CHKERRQ(ierr); 2297 PetscStackPop; 2298 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 2299 ierr = MatDestroy(&vecMat);CHKERRQ(ierr); 2300 ierr = ISRestoreIndices(patch->allCells, &cellsArray);CHKERRQ(ierr); 2301 ierr = ISRestoreIndices(dofMap, &dofMapArray);CHKERRQ(ierr); 2302 ierr = ISDestroy(&dofMap);CHKERRQ(ierr); 2303 2304 if (patch->usercomputeopintfacet) { 2305 PetscInt nIntFacets; 2306 IS intFacetsIS; 2307 const PetscInt *intFacetsArray = NULL; 2308 if (!patch->allIntFacets) { 2309 PetscHSetI facets; 2310 PetscHashIter hi; 2311 PetscInt pStart, pEnd, fStart, fEnd; 2312 PetscInt *allIntFacets = NULL; 2313 ierr = PetscHSetICreate(&facets);CHKERRQ(ierr); 2314 ierr = ISGetIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 2315 ierr = PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);CHKERRQ(ierr); 2316 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2317 for (i = pStart; i < pEnd; i++) { 2318 ierr = PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);CHKERRQ(ierr); 2319 ierr = PetscSectionGetOffset(patch->intFacetCounts, i, &offset);CHKERRQ(ierr); 2320 if (nIntFacets <= 0) continue; 2321 for (j = 0; j < nIntFacets; j++) { 2322 PetscHSetIAdd(facets, intFacetsArray[offset + j]);CHKERRQ(ierr); 2323 } 2324 } 2325 ierr = ISRestoreIndices(patch->intFacets, &intFacetsArray);CHKERRQ(ierr); 2326 ierr = PetscHSetIGetSize(facets, &nIntFacets);CHKERRQ(ierr); 2327 ierr = PetscMalloc1(nIntFacets, &allIntFacets);CHKERRQ(ierr); 2328 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 2329 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2330 ierr = PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);CHKERRQ(ierr); 2331 i = 0; 2332 PetscHashIterBegin(facets, hi); 2333 while (!PetscHashIterAtEnd(facets, hi)) { 2334 PetscHashIterGetKey(facets, hi, allIntFacets[i]); 2335 patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i; 2336 PetscHashIterNext(facets, hi); 2337 i++; 2338 } 2339 ierr = PetscHSetIDestroy(&facets);CHKERRQ(ierr); 2340 ierr = ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);CHKERRQ(ierr); 2341 } 2342 ierr = ISGetSize(patch->allIntFacets, &nIntFacets);CHKERRQ(ierr); 2343 if (!patch->intFacetMats) { 2344 ierr = VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);CHKERRQ(ierr); 2345 ierr = VecSetBlockSize(patch->intFacetMats, ndof*2);CHKERRQ(ierr); 2346 } 2347 ierr = VecSet(patch->intFacetMats, 0);CHKERRQ(ierr); 2348 2349 ierr = MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, 2350 (void*)patch->intFacetMats, &vecMat); 2351 ierr = MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);CHKERRQ(ierr); 2352 ierr = ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);CHKERRQ(ierr); 2353 ierr = ISGetIndices(dofMap, &dofMapArray);CHKERRQ(ierr); 2354 ierr = ISGetIndices(patch->allIntFacets, &intFacetsArray);CHKERRQ(ierr); 2355 ierr = ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);CHKERRQ(ierr); 2356 PetscStackPush("PCPatch user callback (interior facets)"); 2357 /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */ 2358 ierr = patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);CHKERRQ(ierr); 2359 PetscStackPop; 2360 ierr = ISDestroy(&intFacetsIS);CHKERRQ(ierr); 2361 ierr = MatDestroy(&vecMat);CHKERRQ(ierr); 2362 ierr = ISRestoreIndices(patch->allIntFacets, &intFacetsArray);CHKERRQ(ierr); 2363 ierr = ISRestoreIndices(dofMap, &dofMapArray);CHKERRQ(ierr); 2364 ierr = ISDestroy(&dofMap);CHKERRQ(ierr); 2365 } 2366 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 2367 2368 PetscFunctionReturn(0); 2369 } 2370 2371 PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype) 2372 { 2373 PC_PATCH *patch = (PC_PATCH *) pc->data; 2374 const PetscScalar *xArray = NULL; 2375 PetscScalar *yArray = NULL; 2376 const PetscInt *gtolArray = NULL; 2377 PetscInt dof, offset, lidx; 2378 PetscErrorCode ierr; 2379 2380 PetscFunctionBeginHot; 2381 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 2382 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 2383 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 2384 if (scattertype == SCATTER_WITHARTIFICIAL) { 2385 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 2386 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);CHKERRQ(ierr); 2387 ierr = ISGetIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 2388 } else if (scattertype == SCATTER_WITHALL) { 2389 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);CHKERRQ(ierr); 2390 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);CHKERRQ(ierr); 2391 ierr = ISGetIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 2392 } else { 2393 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 2394 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 2395 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2396 } 2397 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 2398 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 2399 for (lidx = 0; lidx < dof; ++lidx) { 2400 const PetscInt gidx = gtolArray[offset+lidx]; 2401 2402 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 2403 else yArray[gidx] += xArray[lidx]; /* Reverse */ 2404 } 2405 if (scattertype == SCATTER_WITHARTIFICIAL) { 2406 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 2407 } else if (scattertype == SCATTER_WITHALL) { 2408 ierr = ISRestoreIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 2409 } else { 2410 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2411 } 2412 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 2413 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 2414 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 2415 PetscFunctionReturn(0); 2416 } 2417 2418 static PetscErrorCode PCSetUp_PATCH_Linear(PC pc) 2419 { 2420 PC_PATCH *patch = (PC_PATCH *) pc->data; 2421 const char *prefix; 2422 PetscInt i; 2423 PetscErrorCode ierr; 2424 2425 PetscFunctionBegin; 2426 if (!pc->setupcalled) { 2427 ierr = PetscMalloc1(patch->npatch, &patch->solver);CHKERRQ(ierr); 2428 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 2429 for (i = 0; i < patch->npatch; ++i) { 2430 KSP ksp; 2431 PC subpc; 2432 2433 ierr = KSPCreate(PETSC_COMM_SELF, &ksp);CHKERRQ(ierr); 2434 ierr = KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);CHKERRQ(ierr); 2435 ierr = KSPSetOptionsPrefix(ksp, prefix);CHKERRQ(ierr); 2436 ierr = KSPAppendOptionsPrefix(ksp, "sub_");CHKERRQ(ierr); 2437 ierr = PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);CHKERRQ(ierr); 2438 ierr = KSPGetPC(ksp, &subpc);CHKERRQ(ierr); 2439 ierr = PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);CHKERRQ(ierr); 2440 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);CHKERRQ(ierr); 2441 patch->solver[i] = (PetscObject) ksp; 2442 } 2443 } 2444 if (patch->save_operators) { 2445 if (patch->precomputeElementTensors) { 2446 ierr = PCPatchPrecomputePatchTensors_Private(pc);CHKERRQ(ierr); 2447 } 2448 for (i = 0; i < patch->npatch; ++i) { 2449 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 2450 ierr = PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);CHKERRQ(ierr); 2451 ierr = KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 2452 } 2453 } 2454 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2455 for (i = 0; i < patch->npatch; ++i) { 2456 /* Instead of padding patch->patchUpdate with zeros to get */ 2457 /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */ 2458 /* just get rid of the columns that correspond to the dofs with */ 2459 /* artificial bcs. That's of course fairly inefficient, hopefully we */ 2460 /* can just assemble the rectangular matrix in the first place. */ 2461 Mat matSquare; 2462 IS rowis; 2463 PetscInt dof; 2464 2465 ierr = MatGetSize(patch->mat[i], &dof, NULL);CHKERRQ(ierr); 2466 if (dof == 0) { 2467 patch->matWithArtificial[i] = NULL; 2468 continue; 2469 } 2470 2471 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 2472 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 2473 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 2474 2475 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 2476 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 2477 if(pc->setupcalled) { 2478 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 2479 } else { 2480 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 2481 } 2482 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 2483 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 2484 } 2485 } 2486 PetscFunctionReturn(0); 2487 } 2488 2489 static PetscErrorCode PCSetUp_PATCH(PC pc) 2490 { 2491 PC_PATCH *patch = (PC_PATCH *) pc->data; 2492 PetscInt i; 2493 PetscBool isNonlinear; 2494 PetscErrorCode ierr; 2495 2496 PetscFunctionBegin; 2497 if (!pc->setupcalled) { 2498 PetscInt pStart, pEnd, p; 2499 PetscInt localSize; 2500 2501 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 2502 2503 isNonlinear = patch->isNonlinear; 2504 if (!patch->nsubspaces) { 2505 DM dm; 2506 PetscSection s; 2507 PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 2508 2509 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 2510 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 2511 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 2512 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 2513 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 2514 for (p = pStart; p < pEnd; ++p) { 2515 PetscInt cdof; 2516 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 2517 numGlobalBcs += cdof; 2518 } 2519 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2520 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 2521 for (f = 0; f < Nf; ++f) { 2522 PetscFE fe; 2523 PetscDualSpace sp; 2524 PetscInt cdoff = 0; 2525 2526 ierr = DMGetField(dm, f, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 2527 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 2528 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 2529 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 2530 totNb += Nb[f]; 2531 2532 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 2533 for (c = cStart; c < cEnd; ++c) { 2534 PetscInt *closure = NULL; 2535 PetscInt clSize = 0, cl; 2536 2537 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2538 for (cl = 0; cl < clSize*2; cl += 2) { 2539 const PetscInt p = closure[cl]; 2540 PetscInt fdof, d, foff; 2541 2542 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 2543 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 2544 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 2545 } 2546 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2547 } 2548 if (cdoff != (cEnd-cStart)*Nb[f]) SETERRQ4(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %D for field %D should be Nc (%D) * cellDof (%D)", cdoff, f, cEnd-cStart, Nb[f]); 2549 } 2550 numGlobalBcs = 0; 2551 for (p = pStart; p < pEnd; ++p) { 2552 const PetscInt *ind; 2553 PetscInt off, cdof, d; 2554 2555 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 2556 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 2557 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 2558 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 2559 } 2560 2561 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 2562 for (f = 0; f < Nf; ++f) { 2563 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 2564 } 2565 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 2566 ierr = PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);CHKERRQ(ierr); 2567 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 2568 } 2569 2570 localSize = patch->subspaceOffsets[patch->nsubspaces]; 2571 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);CHKERRQ(ierr); 2572 ierr = VecSetUp(patch->localRHS);CHKERRQ(ierr); 2573 ierr = VecDuplicate(patch->localRHS, &patch->localUpdate);CHKERRQ(ierr); 2574 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 2575 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 2576 2577 /* OK, now build the work vectors */ 2578 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 2579 ierr = PetscMalloc1(patch->npatch, &patch->patchRHS);CHKERRQ(ierr); 2580 ierr = PetscMalloc1(patch->npatch, &patch->patchUpdate);CHKERRQ(ierr); 2581 2582 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2583 ierr = PetscMalloc1(patch->npatch, &patch->patchRHSWithArtificial);CHKERRQ(ierr); 2584 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 2585 } 2586 if (isNonlinear) { 2587 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 2588 } 2589 for (p = pStart; p < pEnd; ++p) { 2590 PetscInt dof; 2591 2592 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 2593 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHS[p-pStart]);CHKERRQ(ierr); 2594 ierr = VecSetUp(patch->patchRHS[p-pStart]);CHKERRQ(ierr); 2595 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 2596 ierr = VecSetUp(patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 2597 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2598 const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL; 2599 PetscInt numPatchDofs, offset; 2600 PetscInt numPatchDofsWithArtificial, offsetWithArtificial; 2601 PetscInt dofWithoutArtificialCounter = 0; 2602 PetscInt *patchWithoutArtificialToWithArtificialArray; 2603 2604 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 2605 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 2606 ierr = VecSetUp(patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 2607 2608 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 2609 /* the index in the patch with all dofs */ 2610 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2611 2612 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 2613 if (numPatchDofs == 0) { 2614 patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL; 2615 continue; 2616 } 2617 2618 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 2619 ierr = ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 2620 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);CHKERRQ(ierr); 2621 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);CHKERRQ(ierr); 2622 2623 ierr = PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);CHKERRQ(ierr); 2624 for (i=0; i<numPatchDofsWithArtificial; i++) { 2625 if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) { 2626 patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i; 2627 dofWithoutArtificialCounter++; 2628 if (dofWithoutArtificialCounter == numPatchDofs) 2629 break; 2630 } 2631 } 2632 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);CHKERRQ(ierr); 2633 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2634 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 2635 } 2636 if (isNonlinear) { 2637 const PetscInt *gtolArray, *gtolArrayWithAll = NULL; 2638 PetscInt numPatchDofs, offset; 2639 PetscInt numPatchDofsWithAll, offsetWithAll; 2640 PetscInt dofWithoutAllCounter = 0; 2641 PetscInt *patchWithoutAllToWithAllArray; 2642 2643 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 2644 /* the index in the patch with all dofs */ 2645 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2646 2647 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 2648 if (numPatchDofs == 0) { 2649 patch->dofMappingWithoutToWithAll[p-pStart] = NULL; 2650 continue; 2651 } 2652 2653 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 2654 ierr = ISGetIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 2655 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);CHKERRQ(ierr); 2656 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);CHKERRQ(ierr); 2657 2658 ierr = PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);CHKERRQ(ierr); 2659 2660 for (i=0; i<numPatchDofsWithAll; i++) { 2661 if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) { 2662 patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i; 2663 dofWithoutAllCounter++; 2664 if (dofWithoutAllCounter == numPatchDofs) 2665 break; 2666 } 2667 } 2668 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);CHKERRQ(ierr); 2669 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2670 ierr = ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 2671 } 2672 } 2673 if (patch->save_operators) { 2674 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 2675 for (i = 0; i < patch->npatch; ++i) { 2676 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);CHKERRQ(ierr); 2677 } 2678 } 2679 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 2680 2681 /* If desired, calculate weights for dof multiplicity */ 2682 if (patch->partition_of_unity) { 2683 PetscScalar *input = NULL; 2684 PetscScalar *output = NULL; 2685 Vec global; 2686 2687 ierr = VecDuplicate(patch->localRHS, &patch->dof_weights);CHKERRQ(ierr); 2688 if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) { 2689 for (i = 0; i < patch->npatch; ++i) { 2690 PetscInt dof; 2691 2692 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 2693 if (dof <= 0) continue; 2694 ierr = VecSet(patch->patchRHS[i], 1.0);CHKERRQ(ierr); 2695 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 2696 } 2697 } else { 2698 /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */ 2699 ierr = VecSet(patch->dof_weights, 1.0);CHKERRQ(ierr); 2700 } 2701 2702 VecDuplicate(patch->dof_weights, &global); 2703 VecSet(global, 0.); 2704 2705 ierr = VecGetArray(patch->dof_weights, &input);CHKERRQ(ierr); 2706 ierr = VecGetArray(global, &output);CHKERRQ(ierr); 2707 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2708 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2709 ierr = VecRestoreArray(patch->dof_weights, &input);CHKERRQ(ierr); 2710 ierr = VecRestoreArray(global, &output);CHKERRQ(ierr); 2711 2712 ierr = VecReciprocal(global);CHKERRQ(ierr); 2713 2714 ierr = VecGetArray(patch->dof_weights, &output);CHKERRQ(ierr); 2715 ierr = VecGetArray(global, &input);CHKERRQ(ierr); 2716 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2717 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2718 ierr = VecRestoreArray(patch->dof_weights, &output);CHKERRQ(ierr); 2719 ierr = VecRestoreArray(global, &input);CHKERRQ(ierr); 2720 ierr = VecDestroy(&global);CHKERRQ(ierr); 2721 } 2722 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) { 2723 ierr = PetscMalloc1(patch->npatch, &patch->matWithArtificial);CHKERRQ(ierr); 2724 } 2725 } 2726 ierr = (*patch->setupsolver)(pc);CHKERRQ(ierr); 2727 PetscFunctionReturn(0); 2728 } 2729 2730 static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y) 2731 { 2732 PC_PATCH *patch = (PC_PATCH *) pc->data; 2733 KSP ksp = (KSP) patch->solver[i]; 2734 PetscErrorCode ierr; 2735 2736 PetscFunctionBegin; 2737 if (!patch->save_operators) { 2738 Mat mat; 2739 2740 ierr = PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);CHKERRQ(ierr); 2741 /* Populate operator here. */ 2742 ierr = PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);CHKERRQ(ierr); 2743 ierr = KSPSetOperators(ksp, mat, mat);CHKERRQ(ierr); 2744 /* Drop reference so the KSPSetOperators below will blow it away. */ 2745 ierr = MatDestroy(&mat);CHKERRQ(ierr); 2746 } 2747 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2748 if (!ksp->setfromoptionscalled) { 2749 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 2750 } 2751 ierr = KSPSolve(ksp, x, y);CHKERRQ(ierr); 2752 ierr = KSPCheckSolve(ksp, pc, y);CHKERRQ(ierr); 2753 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2754 if (!patch->save_operators) { 2755 PC pc; 2756 ierr = KSPSetOperators(ksp, NULL, NULL);CHKERRQ(ierr); 2757 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 2758 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 2759 ierr = PCReset(pc);CHKERRQ(ierr); 2760 } 2761 PetscFunctionReturn(0); 2762 } 2763 2764 static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart) 2765 { 2766 PC_PATCH *patch = (PC_PATCH *) pc->data; 2767 Mat multMat; 2768 PetscErrorCode ierr; 2769 2770 PetscFunctionBegin; 2771 2772 if (patch->save_operators) { 2773 multMat = patch->matWithArtificial[i]; 2774 } else { 2775 /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/ 2776 Mat matSquare; 2777 PetscInt dof; 2778 IS rowis; 2779 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 2780 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 2781 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 2782 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 2783 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 2784 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat); CHKERRQ(ierr); 2785 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 2786 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 2787 } 2788 ierr = MatMult(multMat, patch->patchUpdate[i], patch->patchRHSWithArtificial[i]); CHKERRQ(ierr); 2789 ierr = VecScale(patch->patchRHSWithArtificial[i], -1.0); CHKERRQ(ierr); 2790 ierr = PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial[i], patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL); CHKERRQ(ierr); 2791 if (!patch->save_operators) { 2792 ierr = MatDestroy(&multMat); CHKERRQ(ierr); 2793 } 2794 PetscFunctionReturn(0); 2795 } 2796 2797 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 2798 { 2799 PC_PATCH *patch = (PC_PATCH *) pc->data; 2800 const PetscScalar *globalRHS = NULL; 2801 PetscScalar *localRHS = NULL; 2802 PetscScalar *globalUpdate = NULL; 2803 const PetscInt *bcNodes = NULL; 2804 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 2805 PetscInt start[2] = {0, 0}; 2806 PetscInt end[2] = {-1, -1}; 2807 const PetscInt inc[2] = {1, -1}; 2808 const PetscScalar *localUpdate; 2809 const PetscInt *iterationSet; 2810 PetscInt pStart, numBcs, n, sweep, bc, j; 2811 PetscErrorCode ierr; 2812 2813 PetscFunctionBegin; 2814 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2815 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 2816 /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */ 2817 end[0] = patch->npatch; 2818 start[1] = patch->npatch-1; 2819 if (patch->user_patches) { 2820 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 2821 start[1] = end[0] - 1; 2822 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 2823 } 2824 /* Scatter from global space into overlapped local spaces */ 2825 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2826 ierr = VecGetArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2827 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2828 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2829 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2830 ierr = VecRestoreArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2831 2832 ierr = VecSet(patch->localUpdate, 0.0);CHKERRQ(ierr); 2833 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 2834 for (sweep = 0; sweep < nsweep; sweep++) { 2835 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 2836 PetscInt i = patch->user_patches ? iterationSet[j] : j; 2837 PetscInt start, len; 2838 2839 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 2840 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 2841 /* TODO: Squash out these guys in the setup as well. */ 2842 if (len <= 0) continue; 2843 /* TODO: Do we need different scatters for X and Y? */ 2844 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS[i], INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);CHKERRQ(ierr); 2845 ierr = (*patch->applysolver)(pc, i, patch->patchRHS[i], patch->patchUpdate[i]);CHKERRQ(ierr); 2846 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate[i], patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 2847 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2848 ierr = (*patch->updatemultiplicative)(pc, i, pStart);CHKERRQ(ierr); 2849 } 2850 } 2851 } 2852 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 2853 /* XXX: should we do this on the global vector? */ 2854 if (patch->partition_of_unity) { 2855 ierr = VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);CHKERRQ(ierr); 2856 } 2857 /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */ 2858 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 2859 ierr = VecGetArray(y, &globalUpdate);CHKERRQ(ierr); 2860 ierr = VecGetArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2861 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2862 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2863 ierr = VecRestoreArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2864 2865 /* Now we need to send the global BC values through */ 2866 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2867 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 2868 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2869 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 2870 for (bc = 0; bc < numBcs; ++bc) { 2871 const PetscInt idx = bcNodes[bc]; 2872 if (idx < n) globalUpdate[idx] = globalRHS[idx]; 2873 } 2874 2875 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2876 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2877 ierr = VecRestoreArray(y, &globalUpdate);CHKERRQ(ierr); 2878 2879 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 2880 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2881 PetscFunctionReturn(0); 2882 } 2883 2884 static PetscErrorCode PCReset_PATCH_Linear(PC pc) 2885 { 2886 PC_PATCH *patch = (PC_PATCH *) pc->data; 2887 PetscInt i; 2888 PetscErrorCode ierr; 2889 2890 PetscFunctionBegin; 2891 if (patch->solver) { 2892 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset((KSP) patch->solver[i]);CHKERRQ(ierr);} 2893 } 2894 PetscFunctionReturn(0); 2895 } 2896 2897 static PetscErrorCode PCReset_PATCH(PC pc) 2898 { 2899 PC_PATCH *patch = (PC_PATCH *) pc->data; 2900 PetscInt i; 2901 PetscErrorCode ierr; 2902 2903 PetscFunctionBegin; 2904 2905 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 2906 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 2907 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 2908 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 2909 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 2910 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 2911 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 2912 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 2913 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 2914 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 2915 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 2916 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 2917 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 2918 ierr = PetscSectionDestroy(&patch->gtolCountsWithArtificial);CHKERRQ(ierr); 2919 ierr = ISDestroy(&patch->gtolWithArtificial);CHKERRQ(ierr); 2920 ierr = ISDestroy(&patch->dofsWithArtificial);CHKERRQ(ierr); 2921 ierr = ISDestroy(&patch->offsWithArtificial);CHKERRQ(ierr); 2922 ierr = PetscSectionDestroy(&patch->gtolCountsWithAll);CHKERRQ(ierr); 2923 ierr = ISDestroy(&patch->gtolWithAll);CHKERRQ(ierr); 2924 ierr = ISDestroy(&patch->dofsWithAll);CHKERRQ(ierr); 2925 ierr = ISDestroy(&patch->offsWithAll);CHKERRQ(ierr); 2926 ierr = VecDestroy(&patch->cellMats);CHKERRQ(ierr); 2927 ierr = VecDestroy(&patch->intFacetMats);CHKERRQ(ierr); 2928 ierr = ISDestroy(&patch->allCells);CHKERRQ(ierr); 2929 2930 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 2931 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 2932 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 2933 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 2934 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 2935 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 2936 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 2937 2938 ierr = (*patch->resetsolver)(pc);CHKERRQ(ierr); 2939 2940 if (patch->subspaces_to_exclude) { 2941 PetscHSetIDestroy(&patch->subspaces_to_exclude); 2942 } 2943 2944 ierr = VecDestroy(&patch->localRHS);CHKERRQ(ierr); 2945 ierr = VecDestroy(&patch->localUpdate);CHKERRQ(ierr); 2946 if (patch->patchRHS) { 2947 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHS[i]);CHKERRQ(ierr);} 2948 ierr = PetscFree(patch->patchRHS);CHKERRQ(ierr); 2949 } 2950 if (patch->patchUpdate) { 2951 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchUpdate[i]);CHKERRQ(ierr);} 2952 ierr = PetscFree(patch->patchUpdate);CHKERRQ(ierr); 2953 } 2954 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 2955 if (patch->patch_dof_weights) { 2956 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 2957 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 2958 } 2959 if (patch->mat) { 2960 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 2961 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 2962 } 2963 if (patch->matWithArtificial) { 2964 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->matWithArtificial[i]);CHKERRQ(ierr);} 2965 ierr = PetscFree(patch->matWithArtificial);CHKERRQ(ierr); 2966 } 2967 if (patch->patchRHSWithArtificial) { 2968 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHSWithArtificial[i]);CHKERRQ(ierr);} 2969 ierr = PetscFree(patch->patchRHSWithArtificial);CHKERRQ(ierr); 2970 } 2971 if(patch->dofMappingWithoutToWithArtificial) { 2972 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);CHKERRQ(ierr);} 2973 ierr = PetscFree(patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 2974 2975 } 2976 if(patch->dofMappingWithoutToWithAll) { 2977 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithAll[i]);CHKERRQ(ierr);} 2978 ierr = PetscFree(patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 2979 2980 } 2981 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 2982 if (patch->userIS) { 2983 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 2984 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 2985 } 2986 ierr = PetscFree(patch->precomputedTensorLocations);CHKERRQ(ierr); 2987 ierr = PetscFree(patch->precomputedIntFacetTensorLocations);CHKERRQ(ierr); 2988 2989 patch->bs = 0; 2990 patch->cellNodeMap = NULL; 2991 patch->nsubspaces = 0; 2992 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 2993 2994 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 2995 PetscFunctionReturn(0); 2996 } 2997 2998 static PetscErrorCode PCDestroy_PATCH_Linear(PC pc) 2999 { 3000 PC_PATCH *patch = (PC_PATCH *) pc->data; 3001 PetscInt i; 3002 PetscErrorCode ierr; 3003 3004 PetscFunctionBegin; 3005 if (patch->solver) { 3006 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy((KSP *) &patch->solver[i]);CHKERRQ(ierr);} 3007 ierr = PetscFree(patch->solver);CHKERRQ(ierr); 3008 } 3009 PetscFunctionReturn(0); 3010 } 3011 3012 static PetscErrorCode PCDestroy_PATCH(PC pc) 3013 { 3014 PC_PATCH *patch = (PC_PATCH *) pc->data; 3015 PetscErrorCode ierr; 3016 3017 PetscFunctionBegin; 3018 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 3019 ierr = (*patch->destroysolver)(pc);CHKERRQ(ierr); 3020 ierr = PetscFree(pc->data);CHKERRQ(ierr); 3021 PetscFunctionReturn(0); 3022 } 3023 3024 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 3025 { 3026 PC_PATCH *patch = (PC_PATCH *) pc->data; 3027 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 3028 char sub_mat_type[PETSC_MAX_PATH_LEN]; 3029 char option[PETSC_MAX_PATH_LEN]; 3030 const char *prefix; 3031 PetscBool flg, dimflg, codimflg; 3032 MPI_Comm comm; 3033 PetscInt *ifields, nfields, k; 3034 PetscErrorCode ierr; 3035 PCCompositeType loctype = PC_COMPOSITE_ADDITIVE; 3036 3037 PetscFunctionBegin; 3038 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 3039 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 3040 ierr = PetscOptionsHead(PetscOptionsObject, "Patch solver options");CHKERRQ(ierr); 3041 3042 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname);CHKERRQ(ierr); 3043 ierr = PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 3044 3045 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname); 3046 ierr = PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg);CHKERRQ(ierr); 3047 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname); 3048 ierr = PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 3049 3050 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);CHKERRQ(ierr); 3051 ierr = PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);CHKERRQ(ierr); 3052 if(flg) { ierr = PCPatchSetLocalComposition(pc, loctype);CHKERRQ(ierr);} 3053 3054 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);CHKERRQ(ierr); 3055 ierr = PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);CHKERRQ(ierr); 3056 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);CHKERRQ(ierr); 3057 ierr = PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);CHKERRQ(ierr); 3058 if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension"); 3059 3060 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname);CHKERRQ(ierr); 3061 ierr = PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 3062 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 3063 3064 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname);CHKERRQ(ierr); 3065 ierr = PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 3066 3067 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname);CHKERRQ(ierr); 3068 ierr = PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 3069 3070 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname);CHKERRQ(ierr); 3071 ierr = PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg);CHKERRQ(ierr); 3072 3073 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname);CHKERRQ(ierr); 3074 ierr = PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);CHKERRQ(ierr); 3075 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 3076 3077 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname);CHKERRQ(ierr); 3078 ierr = PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 3079 3080 /* If the user has set the number of subspaces, use that for the buffer size, 3081 otherwise use a large number */ 3082 if (patch->nsubspaces <= 0) { 3083 nfields = 128; 3084 } else { 3085 nfields = patch->nsubspaces; 3086 } 3087 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 3088 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);CHKERRQ(ierr); 3089 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);CHKERRQ(ierr); 3090 if (flg && (patchConstructionType == PC_PATCH_USER)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point"); 3091 if (flg) { 3092 PetscHSetIClear(patch->subspaces_to_exclude); 3093 for (k = 0; k < nfields; k++) { 3094 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 3095 } 3096 } 3097 ierr = PetscFree(ifields);CHKERRQ(ierr); 3098 3099 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);CHKERRQ(ierr); 3100 ierr = PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 3101 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);CHKERRQ(ierr); 3102 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 3103 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);CHKERRQ(ierr); 3104 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);CHKERRQ(ierr); 3105 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);CHKERRQ(ierr); 3106 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);CHKERRQ(ierr); 3107 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);CHKERRQ(ierr); 3108 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 3109 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);CHKERRQ(ierr); 3110 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 3111 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);CHKERRQ(ierr); 3112 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 3113 ierr = PetscOptionsTail();CHKERRQ(ierr); 3114 patch->optionsSet = PETSC_TRUE; 3115 PetscFunctionReturn(0); 3116 } 3117 3118 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 3119 { 3120 PC_PATCH *patch = (PC_PATCH*) pc->data; 3121 KSPConvergedReason reason; 3122 PetscInt i; 3123 PetscErrorCode ierr; 3124 3125 PetscFunctionBegin; 3126 if (!patch->save_operators) { 3127 /* Can't do this here because the sub KSPs don't have an operator attached yet. */ 3128 PetscFunctionReturn(0); 3129 } 3130 for (i = 0; i < patch->npatch; ++i) { 3131 if (!((KSP) patch->solver[i])->setfromoptionscalled) { 3132 ierr = KSPSetFromOptions((KSP) patch->solver[i]);CHKERRQ(ierr); 3133 } 3134 ierr = KSPSetUp((KSP) patch->solver[i]);CHKERRQ(ierr); 3135 ierr = KSPGetConvergedReason((KSP) patch->solver[i], &reason);CHKERRQ(ierr); 3136 if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR; 3137 } 3138 PetscFunctionReturn(0); 3139 } 3140 3141 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 3142 { 3143 PC_PATCH *patch = (PC_PATCH *) pc->data; 3144 PetscViewer sviewer; 3145 PetscBool isascii; 3146 PetscMPIInt rank; 3147 PetscErrorCode ierr; 3148 3149 PetscFunctionBegin; 3150 /* TODO Redo tabbing with set tbas in new style */ 3151 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 3152 if (!isascii) PetscFunctionReturn(0); 3153 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 3154 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 3155 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 3156 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 3157 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");CHKERRQ(ierr); 3158 } else { 3159 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 3160 } 3161 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 3162 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 3163 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 3164 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 3165 if (!patch->precomputeElementTensors) {ierr = PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");CHKERRQ(ierr);} 3166 else {ierr = PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");CHKERRQ(ierr);} 3167 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 3168 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 3169 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 3170 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 3171 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 3172 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 3173 3174 if (patch->isNonlinear) { 3175 ierr = PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");CHKERRQ(ierr); 3176 } else { 3177 ierr = PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");CHKERRQ(ierr); 3178 } 3179 if (patch->solver) { 3180 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 3181 if (!rank) { 3182 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 3183 ierr = PetscObjectView(patch->solver[0], sviewer);CHKERRQ(ierr); 3184 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 3185 } 3186 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 3187 } else { 3188 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 3189 ierr = PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");CHKERRQ(ierr); 3190 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 3191 } 3192 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 3193 PetscFunctionReturn(0); 3194 } 3195 3196 /*MC 3197 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 3198 small block additive preconditioners. Block definition is based on topology from 3199 a DM and equation numbering from a PetscSection. 3200 3201 Options Database Keys: 3202 + -pc_patch_cells_view - Views the process local cell numbers for each patch 3203 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 3204 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 3205 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 3206 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 3207 3208 Level: intermediate 3209 3210 .seealso: PCType, PCCreate(), PCSetType() 3211 M*/ 3212 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 3213 { 3214 PC_PATCH *patch; 3215 PetscErrorCode ierr; 3216 3217 PetscFunctionBegin; 3218 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 3219 3220 if (patch->subspaces_to_exclude) { 3221 PetscHSetIDestroy(&patch->subspaces_to_exclude); 3222 } 3223 PetscHSetICreate(&patch->subspaces_to_exclude); 3224 3225 patch->classname = "pc"; 3226 patch->isNonlinear = PETSC_FALSE; 3227 3228 /* Set some defaults */ 3229 patch->combined = PETSC_FALSE; 3230 patch->save_operators = PETSC_TRUE; 3231 patch->local_composition_type = PC_COMPOSITE_ADDITIVE; 3232 patch->precomputeElementTensors = PETSC_FALSE; 3233 patch->partition_of_unity = PETSC_FALSE; 3234 patch->codim = -1; 3235 patch->dim = -1; 3236 patch->vankadim = -1; 3237 patch->ignoredim = -1; 3238 patch->pardecomp_overlap = 0; 3239 patch->patchconstructop = PCPatchConstruct_Star; 3240 patch->symmetrise_sweep = PETSC_FALSE; 3241 patch->npatch = 0; 3242 patch->userIS = NULL; 3243 patch->optionsSet = PETSC_FALSE; 3244 patch->iterationSet = NULL; 3245 patch->user_patches = PETSC_FALSE; 3246 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 3247 patch->viewPatches = PETSC_FALSE; 3248 patch->viewCells = PETSC_FALSE; 3249 patch->viewPoints = PETSC_FALSE; 3250 patch->viewSection = PETSC_FALSE; 3251 patch->viewMatrix = PETSC_FALSE; 3252 patch->setupsolver = PCSetUp_PATCH_Linear; 3253 patch->applysolver = PCApply_PATCH_Linear; 3254 patch->resetsolver = PCReset_PATCH_Linear; 3255 patch->destroysolver = PCDestroy_PATCH_Linear; 3256 patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear; 3257 patch->dofMappingWithoutToWithArtificial = NULL; 3258 patch->dofMappingWithoutToWithAll = NULL; 3259 3260 pc->data = (void *) patch; 3261 pc->ops->apply = PCApply_PATCH; 3262 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 3263 pc->ops->setup = PCSetUp_PATCH; 3264 pc->ops->reset = PCReset_PATCH; 3265 pc->ops->destroy = PCDestroy_PATCH; 3266 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 3267 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 3268 pc->ops->view = PCView_PATCH; 3269 pc->ops->applyrichardson = 0; 3270 3271 PetscFunctionReturn(0); 3272 } 3273