1 #include <petsc/private/pcpatchimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/dmpleximpl.h> /* For DMPlexComputeJacobian_Patch_Internal() */ 3 #include <petscsf.h> 4 #include <petscbt.h> 5 #include <petscds.h> 6 7 PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Scatter, PC_Patch_Apply, PC_Patch_Prealloc; 8 9 PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format) 10 { 11 PetscErrorCode ierr; 12 13 PetscFunctionBeginHot; 14 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 15 ierr = PetscObjectView(obj, viewer);CHKERRQ(ierr); 16 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 17 PetscFunctionReturn(0); 18 } 19 20 static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHashI ht) 21 { 22 PetscInt starSize; 23 PetscInt *star = NULL, si; 24 PetscErrorCode ierr; 25 26 PetscFunctionBegin; 27 PetscHashIClear(ht); 28 /* To start with, add the point we care about */ 29 PetscHashIAdd(ht, point, 0); 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) {PetscHashIAdd(ht, star[si], 0);} 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, PetscHashI 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 PetscHashIClear(ht); 48 /* To start with, add the point we care about */ 49 PetscHashIAdd(ht, point, 0); 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 PetscHashIAdd(ht, newpoint, 0); 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 /* The user's already set the patches in patch->userIS. Build the hash tables */ 80 static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHashI ht) 81 { 82 PC_PATCH *patch = (PC_PATCH *) vpatch; 83 IS patchis = patch->userIS[point]; 84 PetscInt n; 85 const PetscInt *patchdata; 86 PetscInt pStart, pEnd, i; 87 PetscErrorCode ierr; 88 89 PetscFunctionBegin; 90 PetscHashIClear(ht); 91 ierr = DMPlexGetChart(dm, &pStart, &pEnd); 92 ierr = ISGetLocalSize(patchis, &n);CHKERRQ(ierr); 93 ierr = ISGetIndices(patchis, &patchdata);CHKERRQ(ierr); 94 for (i = 0; i < n; ++i) { 95 const PetscInt ownedpoint = patchdata[i]; 96 97 if (ownedpoint < pStart || ownedpoint >= pEnd) { 98 SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd); 99 } 100 PetscHashIAdd(ht, ownedpoint, 0); 101 } 102 ierr = ISRestoreIndices(patchis, &patchdata);CHKERRQ(ierr); 103 PetscFunctionReturn(0); 104 } 105 106 static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs) 107 { 108 PC_PATCH *patch = (PC_PATCH *) pc->data; 109 PetscInt i; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 if (n == 1 && bs[0] == 1) { 114 patch->defaultSF = sf[0]; 115 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 116 } else { 117 PetscInt allRoots = 0, allLeaves = 0; 118 PetscInt leafOffset = 0; 119 PetscInt *ilocal = NULL; 120 PetscSFNode *iremote = NULL; 121 PetscInt *remoteOffsets = NULL; 122 PetscInt index = 0; 123 PetscHashI rankToIndex; 124 PetscInt numRanks = 0; 125 PetscSFNode *remote = NULL; 126 PetscSF rankSF; 127 PetscInt *ranks = NULL; 128 PetscInt *offsets = NULL; 129 MPI_Datatype contig; 130 PetscHashI ranksUniq; 131 132 /* First figure out how many dofs there are in the concatenated numbering. 133 * allRoots: number of owned global dofs; 134 * allLeaves: number of visible dofs (global + ghosted). 135 */ 136 for (i = 0; i < n; ++i) { 137 PetscInt nroots, nleaves; 138 139 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr); 140 allRoots += nroots * bs[i]; 141 allLeaves += nleaves * bs[i]; 142 } 143 ierr = PetscMalloc1(allLeaves, &ilocal);CHKERRQ(ierr); 144 ierr = PetscMalloc1(allLeaves, &iremote);CHKERRQ(ierr); 145 /* Now build an SF that just contains process connectivity. */ 146 PetscHashICreate(ranksUniq); 147 for (i = 0; i < n; ++i) { 148 const PetscMPIInt *ranks = NULL; 149 PetscInt nranks, j; 150 151 ierr = PetscSFSetUp(sf[i]);CHKERRQ(ierr); 152 ierr = PetscSFGetRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL);CHKERRQ(ierr); 153 /* These are all the ranks who communicate with me. */ 154 for (j = 0; j < nranks; ++j) { 155 PetscHashIAdd(ranksUniq, (PetscInt) ranks[j], 0); 156 } 157 } 158 PetscHashISize(ranksUniq, numRanks); 159 ierr = PetscMalloc1(numRanks, &remote);CHKERRQ(ierr); 160 ierr = PetscMalloc1(numRanks, &ranks);CHKERRQ(ierr); 161 PetscHashIGetKeys(ranksUniq, &index, ranks); 162 163 PetscHashICreate(rankToIndex); 164 for (i = 0; i < numRanks; ++i) { 165 remote[i].rank = ranks[i]; 166 remote[i].index = 0; 167 PetscHashIAdd(rankToIndex, ranks[i], i); 168 } 169 ierr = PetscFree(ranks);CHKERRQ(ierr); 170 PetscHashIDestroy(ranksUniq); 171 ierr = PetscSFCreate(PetscObjectComm((PetscObject) pc), &rankSF);CHKERRQ(ierr); 172 ierr = PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 173 ierr = PetscSFSetUp(rankSF);CHKERRQ(ierr); 174 /* OK, use it to communicate the root offset on the remote 175 * processes for each subspace. */ 176 ierr = PetscMalloc1(n, &offsets);CHKERRQ(ierr); 177 ierr = PetscMalloc1(n*numRanks, &remoteOffsets);CHKERRQ(ierr); 178 179 offsets[0] = 0; 180 for (i = 1; i < n; ++i) { 181 PetscInt nroots; 182 183 ierr = PetscSFGetGraph(sf[i-1], &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 184 offsets[i] = offsets[i-1] + nroots*bs[i-1]; 185 } 186 /* Offsets are the offsets on the current process of the 187 * global dof numbering for the subspaces. */ 188 ierr = MPI_Type_contiguous(n, MPIU_INT, &contig);CHKERRQ(ierr); 189 ierr = MPI_Type_commit(&contig);CHKERRQ(ierr); 190 191 ierr = PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 192 ierr = PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 193 ierr = MPI_Type_free(&contig);CHKERRQ(ierr); 194 ierr = PetscFree(offsets);CHKERRQ(ierr); 195 ierr = PetscSFDestroy(&rankSF);CHKERRQ(ierr); 196 /* Now remoteOffsets contains the offsets on the remote 197 * processes who communicate with me. So now we can 198 * concatenate the list of SFs into a single one. */ 199 index = 0; 200 for (i = 0; i < n; ++i) { 201 const PetscSFNode *remote = NULL; 202 const PetscInt *local = NULL; 203 PetscInt nroots, nleaves, j; 204 205 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);CHKERRQ(ierr); 206 for (j = 0; j < nleaves; ++j) { 207 PetscInt rank = remote[j].rank; 208 PetscInt idx, rootOffset, k; 209 210 PetscHashIMap(rankToIndex, rank, idx); 211 if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?"); 212 /* Offset on given rank for ith subspace */ 213 rootOffset = remoteOffsets[n*idx + i]; 214 for (k = 0; k < bs[i]; ++k) { 215 ilocal[index] = local[j]*bs[i] + k + leafOffset; 216 iremote[index].rank = remote[j].rank; 217 iremote[index].index = remote[j].index*bs[i] + k + rootOffset; 218 ++index; 219 } 220 } 221 leafOffset += nleaves * bs[i]; 222 } 223 PetscHashIDestroy(rankToIndex); 224 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 225 ierr = PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->defaultSF);CHKERRQ(ierr); 226 ierr = PetscSFSetGraph(patch->defaultSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);CHKERRQ(ierr); 227 } 228 PetscFunctionReturn(0); 229 } 230 231 /* TODO: Docs */ 232 PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim) 233 { 234 PC_PATCH *patch = (PC_PATCH *) pc->data; 235 PetscFunctionBegin; 236 patch->ignoredim = dim; 237 PetscFunctionReturn(0); 238 } 239 240 /* TODO: Docs */ 241 PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim) 242 { 243 PC_PATCH *patch = (PC_PATCH *) pc->data; 244 PetscFunctionBegin; 245 *dim = patch->ignoredim; 246 PetscFunctionReturn(0); 247 } 248 249 /* TODO: Docs */ 250 PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg) 251 { 252 PC_PATCH *patch = (PC_PATCH *) pc->data; 253 PetscFunctionBegin; 254 patch->save_operators = flg; 255 PetscFunctionReturn(0); 256 } 257 258 /* TODO: Docs */ 259 PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg) 260 { 261 PC_PATCH *patch = (PC_PATCH *) pc->data; 262 PetscFunctionBegin; 263 *flg = patch->save_operators; 264 PetscFunctionReturn(0); 265 } 266 267 /* TODO: Docs */ 268 PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg) 269 { 270 PC_PATCH *patch = (PC_PATCH *) pc->data; 271 PetscFunctionBegin; 272 patch->partition_of_unity = flg; 273 PetscFunctionReturn(0); 274 } 275 276 /* TODO: Docs */ 277 PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg) 278 { 279 PC_PATCH *patch = (PC_PATCH *) pc->data; 280 PetscFunctionBegin; 281 *flg = patch->partition_of_unity; 282 PetscFunctionReturn(0); 283 } 284 285 /* TODO: Docs */ 286 PetscErrorCode PCPatchSetMultiplicative(PC pc, PetscBool flg) 287 { 288 PC_PATCH *patch = (PC_PATCH *) pc->data; 289 PetscFunctionBegin; 290 patch->multiplicative = flg; 291 PetscFunctionReturn(0); 292 } 293 294 /* TODO: Docs */ 295 PetscErrorCode PCPatchGetMultiplicative(PC pc, PetscBool *flg) 296 { 297 PC_PATCH *patch = (PC_PATCH *) pc->data; 298 PetscFunctionBegin; 299 *flg = patch->multiplicative; 300 PetscFunctionReturn(0); 301 } 302 303 /* TODO: Docs */ 304 PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type) 305 { 306 PC_PATCH *patch = (PC_PATCH *) pc->data; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 if (patch->sub_mat_type) {ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr);} 311 ierr = PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 /* TODO: Docs */ 316 PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type) 317 { 318 PC_PATCH *patch = (PC_PATCH *) pc->data; 319 PetscFunctionBegin; 320 *sub_mat_type = patch->sub_mat_type; 321 PetscFunctionReturn(0); 322 } 323 324 /* TODO: Docs */ 325 PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering) 326 { 327 PC_PATCH *patch = (PC_PATCH *) pc->data; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 patch->cellNumbering = cellNumbering; 332 ierr = PetscObjectReference((PetscObject) cellNumbering);CHKERRQ(ierr); 333 PetscFunctionReturn(0); 334 } 335 336 /* TODO: Docs */ 337 PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering) 338 { 339 PC_PATCH *patch = (PC_PATCH *) pc->data; 340 PetscFunctionBegin; 341 *cellNumbering = patch->cellNumbering; 342 PetscFunctionReturn(0); 343 } 344 345 /* TODO: Docs */ 346 PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 347 { 348 PC_PATCH *patch = (PC_PATCH *) pc->data; 349 350 PetscFunctionBegin; 351 patch->ctype = ctype; 352 switch (ctype) { 353 case PC_PATCH_STAR: 354 patch->patchconstructop = PCPatchConstruct_Star; 355 break; 356 case PC_PATCH_VANKA: 357 patch->patchconstructop = PCPatchConstruct_Vanka; 358 break; 359 case PC_PATCH_USER: 360 case PC_PATCH_PYTHON: 361 patch->user_patches = PETSC_TRUE; 362 patch->patchconstructop = PCPatchConstruct_User; 363 patch->userpatchconstructionop = func; 364 patch->userpatchconstructctx = ctx; 365 break; 366 default: 367 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 368 } 369 PetscFunctionReturn(0); 370 } 371 372 /* TODO: Docs */ 373 PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx) 374 { 375 PC_PATCH *patch = (PC_PATCH *) pc->data; 376 377 PetscFunctionBegin; 378 *ctype = patch->ctype; 379 switch (patch->ctype) { 380 case PC_PATCH_STAR: 381 case PC_PATCH_VANKA: 382 break; 383 case PC_PATCH_USER: 384 case PC_PATCH_PYTHON: 385 *func = patch->userpatchconstructionop; 386 *ctx = patch->userpatchconstructctx; 387 break; 388 default: 389 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 390 } 391 PetscFunctionReturn(0); 392 } 393 394 /* TODO: Docs */ 395 PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 396 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 397 { 398 PC_PATCH *patch = (PC_PATCH *) pc->data; 399 DM dm; 400 PetscSF *sfs; 401 PetscInt cStart, cEnd, i, j; 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 406 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 407 ierr = PetscMalloc1(nsubspaces, &sfs);CHKERRQ(ierr); 408 ierr = PetscMalloc1(nsubspaces, &patch->dofSection);CHKERRQ(ierr); 409 ierr = PetscMalloc1(nsubspaces, &patch->bs);CHKERRQ(ierr); 410 ierr = PetscMalloc1(nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 411 ierr = PetscMalloc1(nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 412 ierr = PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 413 414 patch->nsubspaces = nsubspaces; 415 patch->totalDofsPerCell = 0; 416 for (i = 0; i < nsubspaces; ++i) { 417 ierr = DMGetDefaultSection(dms[i], &patch->dofSection[i]);CHKERRQ(ierr); 418 ierr = PetscObjectReference((PetscObject) patch->dofSection[i]);CHKERRQ(ierr); 419 ierr = DMGetDefaultSF(dms[i], &sfs[i]);CHKERRQ(ierr); 420 patch->bs[i] = bs[i]; 421 patch->nodesPerCell[i] = nodesPerCell[i]; 422 patch->totalDofsPerCell += nodesPerCell[i]*bs[i]; 423 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i]*bs[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 424 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]*bs[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 425 patch->subspaceOffsets[i] = subspaceOffsets[i]; 426 } 427 ierr = PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);CHKERRQ(ierr); 428 ierr = PetscFree(sfs);CHKERRQ(ierr); 429 430 patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces]; 431 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 432 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 436 /* TODO: Docs */ 437 PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 438 { 439 PC_PATCH *patch = (PC_PATCH *) pc->data; 440 PetscInt cStart, cEnd, i, j; 441 PetscErrorCode ierr; 442 443 PetscFunctionBegin; 444 patch->combined = PETSC_TRUE; 445 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 446 ierr = DMGetNumFields(dm, &patch->nsubspaces);CHKERRQ(ierr); 447 ierr = PetscCalloc1(patch->nsubspaces, &patch->dofSection);CHKERRQ(ierr); 448 ierr = PetscMalloc1(patch->nsubspaces, &patch->bs);CHKERRQ(ierr); 449 ierr = PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 450 ierr = PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 451 ierr = PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 452 ierr = DMGetDefaultSection(dm, &patch->dofSection[0]);CHKERRQ(ierr); 453 ierr = PetscObjectReference((PetscObject) patch->dofSection[0]);CHKERRQ(ierr); 454 ierr = PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);CHKERRQ(ierr); 455 patch->totalDofsPerCell = 0; 456 for (i = 0; i < patch->nsubspaces; ++i) { 457 patch->bs[i] = 1; 458 patch->nodesPerCell[i] = nodesPerCell[i]; 459 patch->totalDofsPerCell += nodesPerCell[i]; 460 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 461 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 462 } 463 ierr = DMGetDefaultSF(dm, &patch->defaultSF);CHKERRQ(ierr); 464 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 465 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 466 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 /*@C 471 472 PCPatchSetComputeOperator - Set the callback used to compute patch matrices 473 474 Input Parameters: 475 + pc - The PC 476 . func - The callback 477 - ctx - The user context 478 479 Level: advanced 480 481 Note: 482 The callback has signature: 483 + usercomputeop(pc, mat, ncell, cells, n, u, ctx) 484 + pc - The PC 485 + mat - The patch matrix 486 + ncell - The number of cells to integrate over 487 + cells - An array of the cell numbers 488 + n - The size of g2l 489 + g2l - The global to local dof translation table 490 + ctx - The user context 491 and can assume that the matrix entries have been set to zero before the call. 492 493 .seealso: PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 494 @*/ 495 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Mat, PetscInt, const PetscInt [], PetscInt, const PetscInt *, void *), void *ctx) 496 { 497 PC_PATCH *patch = (PC_PATCH *) pc->data; 498 499 PetscFunctionBegin; 500 patch->usercomputeop = func; 501 patch->usercomputectx = ctx; 502 PetscFunctionReturn(0); 503 } 504 505 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 506 on exit, cht contains all the topological entities we need to compute their residuals. 507 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 508 here we assume a standard FE sparsity pattern.*/ 509 /* TODO: Use DMPlexGetAdjacency() */ 510 /* TODO: Look at temp buffer management for GetClosure() */ 511 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHashI ht, PetscHashI cht) 512 { 513 DM dm; 514 PetscHashIIter hi; 515 PetscInt point; 516 PetscInt *star = NULL, *closure = NULL; 517 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 518 PetscErrorCode ierr; 519 520 PetscFunctionBegin; 521 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 522 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 523 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 524 PetscHashIClear(cht); 525 PetscHashIIterBegin(ht, hi); 526 while (!PetscHashIIterAtEnd(ht, hi)) { 527 528 PetscHashIIterGetKey(ht, hi, point); 529 PetscHashIIterNext(ht, hi); 530 531 /* Loop over all the cells that this point connects to */ 532 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 533 for (si = 0; si < starSize*2; si += 2) { 534 const PetscInt ownedpoint = star[si]; 535 /* TODO Check for point in cht before running through closure again */ 536 /* now loop over all entities in the closure of that cell */ 537 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 538 for (ci = 0; ci < closureSize*2; ci += 2) { 539 const PetscInt seenpoint = closure[ci]; 540 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 541 PetscHashIAdd(cht, seenpoint, 0); 542 } 543 } 544 } 545 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 546 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 547 PetscFunctionReturn(0); 548 } 549 550 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 551 { 552 PetscErrorCode ierr; 553 554 PetscFunctionBegin; 555 if (combined) { 556 if (f < 0) { 557 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 558 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 559 } else { 560 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 561 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 562 } 563 } else { 564 if (f < 0) { 565 PC_PATCH *patch = (PC_PATCH *) pc->data; 566 PetscInt fdof, g; 567 568 if (dof) { 569 *dof = 0; 570 for (g = 0; g < patch->nsubspaces; ++g) { 571 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 572 *dof += fdof; 573 } 574 } 575 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 576 } else { 577 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 578 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 579 } 580 } 581 PetscFunctionReturn(0); 582 } 583 584 /* Given a hash table with a set of topological entities (pts), compute the degrees of 585 freedom in global concatenated numbering on those entities. 586 For Vanka smoothing, this needs to do something special: ignore dofs of the 587 constraint subspace on entities that aren't the base entity we're building the patch 588 around. */ 589 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHashI pts, PetscHashI dofs, PetscInt base, PetscInt exclude_subspace) 590 { 591 PC_PATCH *patch = (PC_PATCH *) pc->data; 592 PetscHashIIter hi; 593 PetscInt ldof, loff; 594 PetscInt k, p; 595 PetscErrorCode ierr; 596 597 PetscFunctionBegin; 598 PetscHashIClear(dofs); 599 for (k = 0; k < patch->nsubspaces; ++k) { 600 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 601 PetscInt bs = patch->bs[k]; 602 PetscInt j, l; 603 604 if (k == exclude_subspace) { 605 /* only get this subspace dofs at the base entity, not any others */ 606 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 607 if (0 == ldof) continue; 608 for (j = loff; j < ldof + loff; ++j) { 609 for (l = 0; l < bs; ++l) { 610 PetscInt dof = bs*j + l + subspaceOffset; 611 PetscHashIAdd(dofs, dof, 0); 612 } 613 } 614 continue; /* skip the other dofs of this subspace */ 615 } 616 617 PetscHashIIterBegin(pts, hi); 618 while (!PetscHashIIterAtEnd(pts, hi)) { 619 PetscHashIIterGetKey(pts, hi, p); 620 PetscHashIIterNext(pts, hi); 621 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 622 if (0 == ldof) continue; 623 for (j = loff; j < ldof + loff; ++j) { 624 for (l = 0; l < bs; ++l) { 625 PetscInt dof = bs*j + l + subspaceOffset; 626 PetscHashIAdd(dofs, dof, 0); 627 } 628 } 629 } 630 } 631 PetscFunctionReturn(0); 632 } 633 634 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 635 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHashI A, PetscHashI B, PetscHashI C) 636 { 637 PetscHashIIter hi; 638 PetscInt key, val; 639 PetscBool flg; 640 641 PetscFunctionBegin; 642 PetscHashIClear(C); 643 PetscHashIIterBegin(B, hi); 644 while (!PetscHashIIterAtEnd(B, hi)) { 645 PetscHashIIterGetKeyVal(B, hi, key, val); 646 PetscHashIIterNext(B, hi); 647 PetscHashIHasKey(A, key, flg); 648 if (!flg) {PetscHashIAdd(C, key, val);} 649 } 650 PetscFunctionReturn(0); 651 } 652 653 /* 654 * PCPatchCreateCellPatches - create patches. 655 * 656 * Input Parameters: 657 * + dm - The DMPlex object defining the mesh 658 * 659 * Output Parameters: 660 * + cellCounts - Section with counts of cells around each vertex 661 * . cells - IS of the cell point indices of cells in each patch 662 * . pointCounts - Section with counts of cells around each vertex 663 * - point - IS of the cell point indices of cells in each patch 664 */ 665 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 666 { 667 PC_PATCH *patch = (PC_PATCH *) pc->data; 668 DMLabel ghost = NULL; 669 DM dm, plex; 670 PetscHashI ht, cht; 671 PetscSection cellCounts, pointCounts; 672 PetscInt *cellsArray, *pointsArray; 673 PetscInt numCells, numPoints; 674 const PetscInt *leaves; 675 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 676 PetscBool isFiredrake; 677 PetscErrorCode ierr; 678 679 PetscFunctionBegin; 680 /* Used to keep track of the cells in the patch. */ 681 PetscHashICreate(ht); 682 PetscHashICreate(cht); 683 684 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 685 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 686 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 687 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 688 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 689 690 if (patch->user_patches) { 691 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 692 vStart = 0; vEnd = patch->npatch; 693 } else if (patch->codim < 0) { 694 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 695 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 696 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 697 patch->npatch = vEnd - vStart; 698 699 /* These labels mark the owned points. We only create patches around points that this process owns. */ 700 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 701 if (isFiredrake) { 702 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 703 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 704 } else { 705 PetscSF sf; 706 707 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 708 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 709 nleaves = PetscMax(nleaves, 0); 710 } 711 712 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 713 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 714 cellCounts = patch->cellCounts; 715 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 716 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 717 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 718 pointCounts = patch->pointCounts; 719 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 720 /* Count cells and points in the patch surrounding each entity */ 721 for (v = vStart; v < vEnd; ++v) { 722 PetscHashIIter hi; 723 PetscInt chtSize, loc = -1; 724 PetscBool flg; 725 726 if (!patch->user_patches) { 727 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 728 else {ierr = PetscFindInt(v, nleaves, leaves, &loc); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 729 /* Not an owned entity, don't make a cell patch. */ 730 if (flg) continue; 731 } 732 733 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 734 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 735 PetscHashISize(cht, chtSize); 736 /* empty patch, continue */ 737 if (chtSize == 0) continue; 738 739 /* safe because size(cht) > 0 from above */ 740 PetscHashIIterBegin(cht, hi); 741 while (!PetscHashIIterAtEnd(cht, hi)) { 742 PetscInt point, pdof; 743 744 PetscHashIIterGetKey(cht, hi, point); 745 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 746 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 747 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 748 PetscHashIIterNext(cht, hi); 749 } 750 } 751 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 752 753 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 754 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 755 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 756 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 757 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 758 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 759 760 /* Now that we know how much space we need, run through again and actually remember the cells. */ 761 for (v = vStart; v < vEnd; v++ ) { 762 PetscHashIIter hi; 763 PetscInt dof, off, cdof, coff, pdof, n = 0, cn = 0; 764 765 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 766 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 767 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 768 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 769 if (dof <= 0) continue; 770 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 771 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 772 PetscHashIIterBegin(cht, hi); 773 while (!PetscHashIIterAtEnd(cht, hi)) { 774 PetscInt point; 775 776 PetscHashIIterGetKey(cht, hi, point); 777 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 778 if (pdof) {pointsArray[off + n++] = point;} 779 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 780 PetscHashIIterNext(cht, hi); 781 } 782 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); 783 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); 784 } 785 PetscHashIDestroy(ht); 786 PetscHashIDestroy(cht); 787 ierr = DMDestroy(&plex);CHKERRQ(ierr); 788 789 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 790 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 791 if (patch->viewCells) { 792 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 793 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 794 } 795 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 796 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 797 if (patch->viewPoints) { 798 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 799 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 800 } 801 PetscFunctionReturn(0); 802 } 803 804 /* 805 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 806 * 807 * Input Parameters: 808 * + dm - The DMPlex object defining the mesh 809 * . cellCounts - Section with counts of cells around each vertex 810 * . cells - IS of the cell point indices of cells in each patch 811 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 812 * . nodesPerCell - number of nodes per cell. 813 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 814 * 815 * Output Parameters: 816 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 817 * . gtolCounts - Section with counts of dofs per cell patch 818 * - gtol - IS mapping from global dofs to local dofs for each patch. 819 */ 820 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 821 { 822 PC_PATCH *patch = (PC_PATCH *) pc->data; 823 PetscSection cellCounts = patch->cellCounts; 824 PetscSection pointCounts = patch->pointCounts; 825 PetscSection gtolCounts; 826 IS cells = patch->cells; 827 IS points = patch->points; 828 PetscSection cellNumbering = patch->cellNumbering; 829 PetscInt Nf = patch->nsubspaces; 830 PetscInt numCells, numPoints; 831 PetscInt numDofs; 832 PetscInt numGlobalDofs; 833 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 834 PetscInt vStart, vEnd, v; 835 const PetscInt *cellsArray, *pointsArray; 836 PetscInt *newCellsArray = NULL; 837 PetscInt *dofsArray = NULL; 838 PetscInt *offsArray = NULL; 839 PetscInt *asmArray = NULL; 840 PetscInt *globalDofsArray = NULL; 841 PetscInt globalIndex = 0; 842 PetscInt key = 0; 843 PetscInt asmKey = 0; 844 PetscHashI ht; 845 PetscInt pStart, pEnd, p; 846 PetscErrorCode ierr; 847 848 PetscFunctionBegin; 849 /* dofcounts section is cellcounts section * dofPerCell */ 850 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 851 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 852 numDofs = numCells * totalDofsPerCell; 853 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 854 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 855 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 856 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 857 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 858 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 859 gtolCounts = patch->gtolCounts; 860 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 861 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 862 863 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 864 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 865 PetscHashICreate(ht); 866 for (v = vStart; v < vEnd; ++v) { 867 PetscInt localIndex = 0; 868 PetscInt dof, off, i, j, k, l; 869 870 PetscHashIClear(ht); 871 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 872 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 873 if (dof <= 0) continue; 874 875 for (k = 0; k < patch->nsubspaces; ++k) { 876 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 877 PetscInt nodesPerCell = patch->nodesPerCell[k]; 878 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 879 PetscInt bs = patch->bs[k]; 880 881 for (i = off; i < off + dof; ++i) { 882 /* Walk over the cells in this patch. */ 883 const PetscInt c = cellsArray[i]; 884 PetscInt cell = c; 885 886 /* TODO Change this to an IS */ 887 if (cellNumbering) { 888 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 889 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 890 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 891 } 892 newCellsArray[i] = cell; 893 for (j = 0; j < nodesPerCell; ++j) { 894 /* For each global dof, map it into contiguous local storage. */ 895 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 896 /* finally, loop over block size */ 897 for (l = 0; l < bs; ++l) { 898 PetscInt localDof; 899 900 PetscHashIMap(ht, globalDof + l, localDof); 901 if (localDof == -1) { 902 localDof = localIndex++; 903 PetscHashIAdd(ht, globalDof + l, localDof); 904 } 905 if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 906 /* And store. */ 907 dofsArray[globalIndex++] = localDof; 908 } 909 } 910 } 911 } 912 /* How many local dofs in this patch? */ 913 PetscHashISize(ht, dof); 914 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 915 } 916 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 917 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 918 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 919 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 920 921 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 922 for (v = vStart; v < vEnd; ++v) { 923 PetscHashIIter hi; 924 PetscInt dof, off, Np, ooff, i, j, k, l; 925 926 PetscHashIClear(ht); 927 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 928 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 929 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 930 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 931 if (dof <= 0) continue; 932 933 for (k = 0; k < patch->nsubspaces; ++k) { 934 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 935 PetscInt nodesPerCell = patch->nodesPerCell[k]; 936 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 937 PetscInt bs = patch->bs[k]; 938 939 for (i = off; i < off + dof; ++i) { 940 /* Reconstruct mapping of global-to-local on this patch. */ 941 const PetscInt c = cellsArray[i]; 942 PetscInt cell = c; 943 944 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 945 for (j = 0; j < nodesPerCell; ++j) { 946 for (l = 0; l < bs; ++l) { 947 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 948 const PetscInt localDof = dofsArray[key]; 949 950 key += 1; 951 PetscHashIAdd(ht, globalDof, localDof); 952 } 953 } 954 } 955 if (dof > 0) { 956 /* Shove it in the output data structure. */ 957 PetscInt goff; 958 959 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 960 PetscHashIIterBegin(ht, hi); 961 while (!PetscHashIIterAtEnd(ht, hi)) { 962 PetscInt globalDof, localDof; 963 964 PetscHashIIterGetKeyVal(ht, hi, globalDof, localDof); 965 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 966 PetscHashIIterNext(ht, hi); 967 } 968 } 969 970 for (p = 0; p < Np; ++p) { 971 const PetscInt point = pointsArray[ooff + p]; 972 PetscInt globalDof, localDof; 973 974 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 975 PetscHashIMap(ht, globalDof, localDof); 976 offsArray[(ooff + p)*Nf + k] = localDof; 977 } 978 } 979 980 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 981 We need to create the dof table laid out cellwise first, then by subspace, 982 as the assembler assembles cell-wise and we need to stuff the different 983 contributions of the different function spaces to the right places. So we loop 984 over cells, then over subspaces. */ 985 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 986 for (i = off; i < off + dof; ++i) { 987 const PetscInt c = cellsArray[i]; 988 PetscInt cell = c; 989 990 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 991 for (k = 0; k < patch->nsubspaces; ++k) { 992 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 993 PetscInt nodesPerCell = patch->nodesPerCell[k]; 994 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 995 PetscInt bs = patch->bs[k]; 996 997 for (j = 0; j < nodesPerCell; ++j) { 998 for (l = 0; l < bs; ++l) { 999 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1000 PetscInt localDof; 1001 1002 PetscHashIMap(ht, globalDof, localDof); 1003 asmArray[asmKey++] = localDof; 1004 } 1005 } 1006 } 1007 } 1008 } 1009 } 1010 if (1 == patch->nsubspaces) {ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr);} 1011 1012 PetscHashIDestroy(ht); 1013 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1014 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1015 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1016 /* Create placeholder section for map from points to patch dofs */ 1017 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1018 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1019 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1020 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1021 for (p = pStart; p < pEnd; ++p) { 1022 PetscInt dof, fdof, f; 1023 1024 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1025 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1026 for (f = 0; f < patch->nsubspaces; ++f) { 1027 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1028 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1029 } 1030 } 1031 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1032 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1033 /* Replace cell indices with firedrake-numbered ones. */ 1034 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1035 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1036 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1037 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1038 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1039 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1040 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1041 PetscFunctionReturn(0); 1042 } 1043 1044 static PetscErrorCode PCPatchCreateCellPatchBCs(PC pc) 1045 { 1046 PC_PATCH *patch = (PC_PATCH *) pc->data; 1047 const PetscInt *bcNodes = NULL; 1048 PetscSection gtolCounts = patch->gtolCounts; 1049 IS gtol = patch->gtol; 1050 DM dm; 1051 PetscInt numBcs; 1052 PetscSection bcCounts; 1053 PetscHashI globalBcs, localBcs, patchDofs; 1054 PetscHashI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 1055 PetscHashIIter hi; 1056 PetscInt *bcsArray = NULL; 1057 PetscInt *multBcsArray = NULL; 1058 const PetscInt *gtolArray; 1059 PetscInt vStart, vEnd, v, i; 1060 PetscErrorCode ierr; 1061 1062 PetscFunctionBegin; 1063 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1064 PetscHashICreate(globalBcs); 1065 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes);CHKERRQ(ierr); 1066 ierr = ISGetSize(patch->ghostBcNodes, &numBcs);CHKERRQ(ierr); 1067 for (i = 0; i < numBcs; ++i) { 1068 /* these are already in concatenated numbering */ 1069 PetscHashIAdd(globalBcs, bcNodes[i], 0); 1070 } 1071 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes);CHKERRQ(ierr); 1072 PetscHashICreate(patchDofs); 1073 PetscHashICreate(localBcs); 1074 PetscHashICreate(ownedpts); 1075 PetscHashICreate(seenpts); 1076 PetscHashICreate(owneddofs); 1077 PetscHashICreate(seendofs); 1078 PetscHashICreate(artificialbcs); 1079 1080 ierr = PetscSectionGetChart(patch->cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 1081 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->bcCounts);CHKERRQ(ierr); 1082 bcCounts = patch->bcCounts; 1083 ierr = PetscSectionSetChart(bcCounts, vStart, vEnd);CHKERRQ(ierr); 1084 ierr = PetscMalloc1(vEnd - vStart, &patch->bcs);CHKERRQ(ierr); 1085 if (patch->multiplicative) {ierr = PetscMalloc1(vEnd - vStart, &patch->multBcs);CHKERRQ(ierr);} 1086 1087 ierr = ISGetIndices(gtol, >olArray);CHKERRQ(ierr); 1088 for (v = vStart; v < vEnd; ++v) { 1089 PetscInt bcIndex = 0; 1090 PetscInt multBcIndex = 0; 1091 PetscInt numBcs, dof, off; 1092 1093 PetscHashIClear(patchDofs); 1094 PetscHashIClear(localBcs); 1095 ierr = PetscSectionGetDof(gtolCounts, v, &dof);CHKERRQ(ierr); 1096 ierr = PetscSectionGetOffset(gtolCounts, v, &off);CHKERRQ(ierr); 1097 1098 if (dof <= 0) { 1099 patch->bcs[v-vStart] = NULL; 1100 if (patch->multiplicative) patch->multBcs[v-vStart] = NULL; 1101 continue; 1102 } 1103 1104 for (i = off; i < off + dof; ++i) { 1105 const PetscInt globalDof = gtolArray[i]; 1106 const PetscInt localDof = i-off; 1107 PetscBool flg; 1108 1109 PetscHashIAdd(patchDofs, globalDof, localDof); 1110 PetscHashIHasKey(globalBcs, globalDof, flg); 1111 if (flg) {PetscHashIAdd(localBcs, localDof, 0);} 1112 } 1113 1114 /* If we're doing multiplicative, make the BC data structures now corresponding solely to actual globally imposed Dirichlet BCs */ 1115 if (patch->multiplicative) { 1116 PetscHashISize(localBcs, numBcs); 1117 ierr = PetscMalloc1(numBcs, &multBcsArray);CHKERRQ(ierr); 1118 PetscHashIGetKeys(localBcs, &multBcIndex, multBcsArray); 1119 ierr = PetscSortInt(numBcs, multBcsArray);CHKERRQ(ierr); 1120 ierr = ISCreateGeneral(PETSC_COMM_SELF, numBcs, multBcsArray, PETSC_OWN_POINTER, &(patch->multBcs[v-vStart]));CHKERRQ(ierr); 1121 } 1122 1123 /* Now figure out the artificial BCs: the set difference of {dofs on entities I see on the patch}\{dofs I am responsible for updating} */ 1124 ierr = patch->patchconstructop((void*) patch, dm, v, ownedpts);CHKERRQ(ierr); 1125 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts);CHKERRQ(ierr); 1126 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, patch->exclude_subspace);CHKERRQ(ierr); 1127 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, -1);CHKERRQ(ierr); 1128 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs);CHKERRQ(ierr); 1129 1130 if (patch->viewPatches) { 1131 PetscHashI globalbcdofs; 1132 MPI_Comm comm; 1133 1134 PetscHashICreate(globalbcdofs); 1135 1136 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 1137 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);CHKERRQ(ierr); 1138 PetscHashIIterBegin(owneddofs, hi); 1139 while (!PetscHashIIterAtEnd(owneddofs, hi)) { 1140 PetscInt globalDof; 1141 1142 PetscHashIIterGetKey(owneddofs, hi, globalDof); 1143 PetscHashIIterNext(owneddofs, hi); 1144 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1145 } 1146 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1147 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);CHKERRQ(ierr); 1148 PetscHashIIterBegin(seendofs, hi); 1149 while (!PetscHashIIterAtEnd(seendofs, hi)) { 1150 PetscInt globalDof; 1151 PetscBool flg; 1152 1153 PetscHashIIterGetKey(seendofs, hi, globalDof); 1154 PetscHashIIterNext(seendofs, hi); 1155 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1156 PetscHashIHasKey(globalBcs, globalDof, flg); 1157 if (flg) {PetscHashIAdd(globalbcdofs, globalDof, 0);} 1158 } 1159 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1160 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v);CHKERRQ(ierr); 1161 PetscHashISize(globalbcdofs, numBcs); 1162 if (numBcs > 0) { 1163 PetscHashIIterBegin(globalbcdofs, hi); 1164 while (!PetscHashIIterAtEnd(globalbcdofs, hi)) { 1165 PetscInt globalDof; 1166 PetscHashIIterGetKey(globalbcdofs, hi, globalDof); 1167 PetscHashIIterNext(globalbcdofs, hi); 1168 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1169 } 1170 } 1171 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1172 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 1173 PetscHashISize(artificialbcs, numBcs); 1174 if (numBcs > 0) { 1175 PetscHashIIterBegin(artificialbcs, hi); 1176 while (!PetscHashIIterAtEnd(artificialbcs, hi)) { 1177 PetscInt globalDof; 1178 PetscHashIIterGetKey(artificialbcs, hi, globalDof); 1179 PetscHashIIterNext(artificialbcs, hi); 1180 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1181 } 1182 } 1183 ierr = PetscSynchronizedPrintf(comm, "\n\n");CHKERRQ(ierr); 1184 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 1185 PetscHashIDestroy(globalbcdofs); 1186 } 1187 1188 PetscHashISize(artificialbcs, numBcs); 1189 if (numBcs > 0) { 1190 PetscHashIIterBegin(artificialbcs, hi); 1191 while (!PetscHashIIterAtEnd(artificialbcs, hi)) { 1192 PetscInt globalDof, localDof; 1193 PetscHashIIterGetKey(artificialbcs, hi, globalDof); 1194 PetscHashIIterNext(artificialbcs, hi); 1195 PetscHashIMap(patchDofs, globalDof, localDof); 1196 if (localDof == -1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Patch %d Didn't find dof %d in patch\n", v, globalDof); 1197 PetscHashIAdd(localBcs, localDof, 0); 1198 } 1199 } 1200 1201 /* OK, now we have a hash table with all the bcs indicated by the artificial and global bcs */ 1202 PetscHashISize(localBcs, numBcs); 1203 ierr = PetscSectionSetDof(bcCounts, v, numBcs);CHKERRQ(ierr); 1204 ierr = PetscMalloc1(numBcs, &bcsArray);CHKERRQ(ierr); 1205 PetscHashIGetKeys(localBcs, &bcIndex, bcsArray); 1206 ierr = PetscSortInt(numBcs, bcsArray);CHKERRQ(ierr); 1207 ierr = ISCreateGeneral(PETSC_COMM_SELF, numBcs, bcsArray, PETSC_OWN_POINTER, &(patch->bcs[v - vStart]));CHKERRQ(ierr); 1208 } 1209 ierr = ISRestoreIndices(gtol, >olArray);CHKERRQ(ierr); 1210 PetscHashIDestroy(artificialbcs); 1211 PetscHashIDestroy(seendofs); 1212 PetscHashIDestroy(owneddofs); 1213 PetscHashIDestroy(seenpts); 1214 PetscHashIDestroy(ownedpts); 1215 PetscHashIDestroy(localBcs); 1216 PetscHashIDestroy(patchDofs); 1217 PetscHashIDestroy(globalBcs); 1218 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 1219 ierr = PetscSectionSetUp(bcCounts);CHKERRQ(ierr); 1220 PetscFunctionReturn(0); 1221 } 1222 1223 static PetscErrorCode PCPatchZeroFillMatrix_Private(Mat mat, const PetscInt ncell, const PetscInt ndof, const PetscInt *dof) 1224 { 1225 const PetscScalar *values = NULL; 1226 PetscInt rows, c, i; 1227 PetscErrorCode ierr; 1228 1229 PetscFunctionBegin; 1230 ierr = PetscCalloc1(ndof*ndof, &values);CHKERRQ(ierr); 1231 for (c = 0; c < ncell; ++c) { 1232 const PetscInt *idx = &dof[ndof*c]; 1233 ierr = MatSetValues(mat, ndof, idx, ndof, idx, values, INSERT_VALUES);CHKERRQ(ierr); 1234 } 1235 ierr = MatGetLocalSize(mat, &rows, NULL);CHKERRQ(ierr); 1236 for (i = 0; i < rows; ++i) { 1237 ierr = MatSetValues(mat, 1, &i, 1, &i, values, INSERT_VALUES);CHKERRQ(ierr); 1238 } 1239 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1240 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1241 ierr = PetscFree(values);CHKERRQ(ierr); 1242 PetscFunctionReturn(0); 1243 } 1244 1245 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat) 1246 { 1247 PC_PATCH *patch = (PC_PATCH *) pc->data; 1248 Vec x, y; 1249 PetscBool flg; 1250 PetscInt csize, rsize; 1251 const char *prefix = NULL; 1252 PetscErrorCode ierr; 1253 1254 PetscFunctionBegin; 1255 x = patch->patchX[point]; 1256 y = patch->patchY[point]; 1257 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1258 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1259 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1260 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1261 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1262 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1263 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1264 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1265 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1266 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1267 /* Sparse patch matrices */ 1268 if (!flg) { 1269 PetscBT bt; 1270 PetscInt *dnnz = NULL; 1271 const PetscInt *dofsArray = NULL; 1272 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1273 1274 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1275 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1276 point += pStart; 1277 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1278 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1279 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1280 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1281 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1282 /* XXX: This uses N^2 bits to store the sparsity pattern on a 1283 * patch. This is probably OK if the patches are not too big, 1284 * but could use quite a bit of memory for planes in 3D. 1285 * Should we switch based on the value of rsize to a 1286 * hash-table (slower, but more memory efficient) approach? */ 1287 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1288 for (c = 0; c < ncell; ++c) { 1289 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1290 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1291 const PetscInt row = idx[i]; 1292 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1293 const PetscInt col = idx[j]; 1294 const PetscInt key = row*rsize + col; 1295 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1296 } 1297 } 1298 } 1299 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1300 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1301 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1302 ierr = PCPatchZeroFillMatrix_Private(*mat, ncell, patch->totalDofsPerCell, &dofsArray[offset*patch->totalDofsPerCell]);CHKERRQ(ierr); 1303 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1304 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1305 } 1306 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Mat J, PetscInt ncell, const PetscInt cells[], PetscInt n, const PetscInt *l2p, void *ctx) 1311 { 1312 PC_PATCH *patch = (PC_PATCH *) pc->data; 1313 DM dm; 1314 PetscSection s; 1315 const PetscInt *parray, *oarray; 1316 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1317 PetscErrorCode ierr; 1318 1319 PetscFunctionBegin; 1320 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1321 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1322 /* Set offset into patch */ 1323 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1324 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1325 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1326 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1327 for (f = 0; f < Nf; ++f) { 1328 for (p = 0; p < Np; ++p) { 1329 const PetscInt point = parray[poff+p]; 1330 PetscInt dof; 1331 1332 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1333 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1334 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1335 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1336 } 1337 } 1338 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1339 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1340 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1341 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 1342 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, 0, ncell, cells, 0.0, 0.0, NULL, NULL, J, J, ctx);CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345 1346 static PetscErrorCode PCPatchComputeOperator_Private(PC pc, Mat mat, Mat multMat, PetscInt point) 1347 { 1348 PC_PATCH *patch = (PC_PATCH *) pc->data; 1349 const PetscInt *dofsArray; 1350 const PetscInt *cellsArray; 1351 PetscInt ncell, offset, pStart, pEnd; 1352 PetscErrorCode ierr; 1353 1354 PetscFunctionBegin; 1355 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1356 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1357 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1358 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1359 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1360 1361 point += pStart; 1362 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1363 1364 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1365 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1366 if (ncell <= 0) { 1367 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1368 PetscFunctionReturn(0); 1369 } 1370 PetscStackPush("PCPatch user callback"); 1371 ierr = patch->usercomputeop(pc, point, mat, ncell, cellsArray + offset, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, patch->usercomputectx);CHKERRQ(ierr); 1372 PetscStackPop; 1373 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1374 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1375 /* TODO Can apply these BCs through the patchSection (would need to changes sizes of patchX/Y etc) */ 1376 /* Apply boundary conditions. Could also do this through the local_to_patch guy. */ 1377 if (patch->multiplicative) { 1378 ierr = MatCopy(mat, multMat, SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1379 ierr = MatZeroRowsColumnsIS(multMat, patch->multBcs[point-pStart], (PetscScalar) 1.0, NULL, NULL);CHKERRQ(ierr); 1380 } 1381 ierr = MatZeroRowsColumnsIS(mat, patch->bcs[point-pStart], (PetscScalar) 1.0, NULL, NULL);CHKERRQ(ierr); 1382 if (patch->viewMatrix) {ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr);} 1383 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1384 PetscFunctionReturn(0); 1385 } 1386 1387 static PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat) 1388 { 1389 PC_PATCH *patch = (PC_PATCH *) pc->data; 1390 const PetscScalar *xArray = NULL; 1391 PetscScalar *yArray = NULL; 1392 const PetscInt *gtolArray = NULL; 1393 PetscInt dof, offset, lidx; 1394 PetscErrorCode ierr; 1395 1396 PetscFunctionBeginHot; 1397 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1398 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 1399 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 1400 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1401 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1402 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1403 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 1404 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 1405 for (lidx = 0; lidx < dof; ++lidx) { 1406 const PetscInt gidx = gtolArray[offset+lidx]; 1407 1408 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 1409 else yArray[gidx] += xArray[lidx]; /* Reverse */ 1410 } 1411 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1412 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 1413 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 1414 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1415 PetscFunctionReturn(0); 1416 } 1417 1418 static PetscErrorCode PCSetUp_PATCH(PC pc) 1419 { 1420 PC_PATCH *patch = (PC_PATCH *) pc->data; 1421 PetscScalar *patchX = NULL; 1422 const PetscInt *bcNodes = NULL; 1423 PetscInt numBcs, i, j; 1424 const char *prefix; 1425 PetscErrorCode ierr; 1426 1427 PetscFunctionBegin; 1428 if (!pc->setupcalled) { 1429 PetscInt pStart, pEnd, p; 1430 PetscInt localSize; 1431 1432 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1433 1434 if (!patch->nsubspaces) { 1435 DM dm; 1436 PetscDS prob; 1437 PetscSection s; 1438 PetscInt cStart, cEnd, c, Nf, f, fsize, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 1439 1440 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1441 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 1442 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1443 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 1444 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 1445 for (p = pStart; p < pEnd; ++p) { 1446 PetscInt cdof; 1447 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1448 numGlobalBcs += cdof; 1449 } 1450 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1451 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1452 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 1453 for (f = 0; f < Nf; ++f) { 1454 PetscFE fe; 1455 PetscDualSpace sp; 1456 PetscSection fs; 1457 PetscInt cdoff = 0; 1458 1459 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1460 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 1461 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1462 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 1463 totNb += Nb[f]; 1464 1465 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr); 1466 ierr = PetscSectionGetStorageSize(fs, &fsize);CHKERRQ(ierr); 1467 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 1468 for (c = cStart; c < cEnd; ++c) { 1469 PetscInt *closure = NULL; 1470 PetscInt clSize = 0, cl; 1471 1472 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1473 for (cl = 0; cl < clSize*2; cl += 2) { 1474 const PetscInt p = closure[cl]; 1475 PetscInt fdof, d, foff; 1476 1477 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 1478 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 1479 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 1480 } 1481 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1482 } 1483 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]); 1484 } 1485 numGlobalBcs = 0; 1486 for (p = pStart; p < pEnd; ++p) { 1487 const PetscInt *ind; 1488 PetscInt off, cdof, d; 1489 1490 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1491 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1492 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 1493 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 1494 } 1495 1496 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 1497 for (f = 0; f < Nf; ++f) { 1498 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 1499 } 1500 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 1501 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 1502 } 1503 1504 localSize = patch->subspaceOffsets[patch->nsubspaces]; 1505 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localX);CHKERRQ(ierr); 1506 ierr = VecSetUp(patch->localX);CHKERRQ(ierr); 1507 ierr = VecDuplicate(patch->localX, &patch->localY);CHKERRQ(ierr); 1508 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 1509 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 1510 ierr = PCPatchCreateCellPatchBCs(pc);CHKERRQ(ierr); 1511 1512 /* OK, now build the work vectors */ 1513 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 1514 ierr = PetscMalloc1(patch->npatch, &patch->patchX);CHKERRQ(ierr); 1515 ierr = PetscMalloc1(patch->npatch, &patch->patchY);CHKERRQ(ierr); 1516 if (patch->partition_of_unity && patch->multiplicative) {ierr = PetscMalloc1(patch->npatch, &patch->patch_dof_weights);CHKERRQ(ierr);} 1517 for (p = pStart; p < pEnd; ++p) { 1518 PetscInt dof; 1519 1520 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1521 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchX[p-pStart]);CHKERRQ(ierr); 1522 ierr = VecSetUp(patch->patchX[p-pStart]);CHKERRQ(ierr); 1523 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchY[p-pStart]);CHKERRQ(ierr); 1524 ierr = VecSetUp(patch->patchY[p-pStart]);CHKERRQ(ierr); 1525 if (patch->partition_of_unity && patch->multiplicative) { 1526 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patch_dof_weights[p-pStart]);CHKERRQ(ierr); 1527 ierr = VecSetUp(patch->patch_dof_weights[p-pStart]);CHKERRQ(ierr); 1528 } 1529 } 1530 ierr = PetscMalloc1(patch->npatch, &patch->ksp);CHKERRQ(ierr); 1531 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1532 for (i = 0; i < patch->npatch; ++i) { 1533 ierr = KSPCreate(PETSC_COMM_SELF, &patch->ksp[i]);CHKERRQ(ierr); 1534 ierr = KSPSetOptionsPrefix(patch->ksp[i], prefix);CHKERRQ(ierr); 1535 ierr = KSPAppendOptionsPrefix(patch->ksp[i], "sub_");CHKERRQ(ierr); 1536 } 1537 if (patch->save_operators) { 1538 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 1539 if (patch->multiplicative) {ierr = PetscMalloc1(patch->npatch, &patch->multMat);CHKERRQ(ierr);} 1540 for (i = 0; i < patch->npatch; ++i) { 1541 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i]);CHKERRQ(ierr); 1542 if (patch->multiplicative) {ierr = MatDuplicate(patch->mat[i], MAT_SHARE_NONZERO_PATTERN, &patch->multMat[i]);CHKERRQ(ierr);} 1543 } 1544 } 1545 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1546 1547 /* If desired, calculate weights for dof multiplicity */ 1548 if (patch->partition_of_unity) { 1549 ierr = VecDuplicate(patch->localX, &patch->dof_weights);CHKERRQ(ierr); 1550 for (i = 0; i < patch->npatch; ++i) { 1551 PetscInt dof; 1552 1553 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 1554 if (dof <= 0) continue; 1555 ierr = VecSet(patch->patchX[i], 1.0);CHKERRQ(ierr); 1556 /* TODO: Do we need different scatters for X and Y? */ 1557 ierr = VecGetArray(patch->patchX[i], &patchX);CHKERRQ(ierr); 1558 /* Apply bcs to patchX (zero entries) */ 1559 ierr = ISGetLocalSize(patch->bcs[i], &numBcs);CHKERRQ(ierr); 1560 ierr = ISGetIndices(patch->bcs[i], &bcNodes);CHKERRQ(ierr); 1561 for (j = 0; j < numBcs; ++j) patchX[bcNodes[j]] = 0; 1562 ierr = ISRestoreIndices(patch->bcs[i], &bcNodes);CHKERRQ(ierr); 1563 ierr = VecRestoreArray(patch->patchX[i], &patchX);CHKERRQ(ierr); 1564 1565 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchX[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1566 } 1567 ierr = VecReciprocal(patch->dof_weights);CHKERRQ(ierr); 1568 if (patch->partition_of_unity && patch->multiplicative) { 1569 for (i = 0; i < patch->npatch; ++i) { 1570 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->dof_weights, patch->patch_dof_weights[i], INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1571 } 1572 } 1573 } 1574 } 1575 if (patch->save_operators) { 1576 for (i = 0; i < patch->npatch; ++i) { 1577 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1578 if (patch->multiplicative) {ierr = PCPatchComputeOperator_Private(pc, patch->mat[i], patch->multMat[i], i);CHKERRQ(ierr);} 1579 else {ierr = PCPatchComputeOperator_Private(pc, patch->mat[i], NULL, i);CHKERRQ(ierr);} 1580 ierr = KSPSetOperators(patch->ksp[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1581 } 1582 } 1583 if (!pc->setupcalled && patch->optionsSet) for (i = 0; i < patch->npatch; ++i) {ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr);} 1584 PetscFunctionReturn(0); 1585 } 1586 1587 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 1588 { 1589 PC_PATCH *patch = (PC_PATCH *) pc->data; 1590 const PetscScalar *globalX = NULL; 1591 PetscScalar *localX = NULL; 1592 PetscScalar *globalY = NULL; 1593 PetscScalar *patchX = NULL; 1594 const PetscInt *bcNodes = NULL; 1595 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 1596 PetscInt start[2] = {0, 0}; 1597 PetscInt end[2] = {-1, -1}; 1598 const PetscInt inc[2] = {1, -1}; 1599 const PetscScalar *localY; 1600 const PetscInt *iterationSet; 1601 PetscInt pStart, numBcs, n, sweep, bc, j; 1602 PetscErrorCode ierr; 1603 1604 PetscFunctionBegin; 1605 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1606 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 1607 end[0] = patch->npatch; 1608 start[1] = patch->npatch-1; 1609 if (patch->user_patches) { 1610 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 1611 start[1] = end[0] - 1; 1612 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 1613 } 1614 /* Scatter from global space into overlapped local spaces */ 1615 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1616 ierr = VecGetArray(patch->localX, &localX);CHKERRQ(ierr); 1617 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1618 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1619 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1620 ierr = VecRestoreArray(patch->localX, &localX);CHKERRQ(ierr); 1621 1622 ierr = VecSet(patch->localY, 0.0);CHKERRQ(ierr); 1623 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 1624 for (sweep = 0; sweep < nsweep; sweep++) { 1625 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 1626 Mat multMat = NULL; 1627 PetscInt i = patch->user_patches ? iterationSet[j] : j; 1628 PetscInt start, len; 1629 1630 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 1631 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 1632 /* TODO: Squash out these guys in the setup as well. */ 1633 if (len <= 0) continue; 1634 /* TODO: Do we need different scatters for X and Y? */ 1635 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localX, patch->patchX[i], INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1636 /* Apply bcs to patchX (zero entries) */ 1637 ierr = VecGetArray(patch->patchX[i], &patchX);CHKERRQ(ierr); 1638 ierr = ISGetLocalSize(patch->bcs[i], &numBcs);CHKERRQ(ierr); 1639 ierr = ISGetIndices(patch->bcs[i], &bcNodes);CHKERRQ(ierr); 1640 for (bc = 0; bc < numBcs; ++bc) patchX[bcNodes[bc]] = 0; 1641 ierr = ISRestoreIndices(patch->bcs[i], &bcNodes);CHKERRQ(ierr); 1642 ierr = VecRestoreArray(patch->patchX[i], &patchX);CHKERRQ(ierr); 1643 if (!patch->save_operators) { 1644 Mat mat; 1645 1646 ierr = PCPatchCreateMatrix_Private(pc, i, &mat);CHKERRQ(ierr); 1647 if (patch->multiplicative) {ierr = MatDuplicate(mat, MAT_SHARE_NONZERO_PATTERN, &multMat);CHKERRQ(ierr);} 1648 /* Populate operator here. */ 1649 ierr = PCPatchComputeOperator_Private(pc, mat, multMat, i);CHKERRQ(ierr); 1650 ierr = KSPSetOperators(patch->ksp[i], mat, mat); 1651 /* Drop reference so the KSPSetOperators below will blow it away. */ 1652 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1653 } else if (patch->multiplicative) { 1654 multMat = patch->multMat[i]; 1655 } 1656 1657 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1658 ierr = KSPSolve(patch->ksp[i], patch->patchX[i], patch->patchY[i]);CHKERRQ(ierr); 1659 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1660 1661 if (!patch->save_operators) { 1662 PC pc; 1663 ierr = KSPSetOperators(patch->ksp[i], NULL, NULL);CHKERRQ(ierr); 1664 ierr = KSPGetPC(patch->ksp[i], &pc);CHKERRQ(ierr); 1665 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 1666 ierr = PCReset(pc);CHKERRQ(ierr); 1667 } 1668 1669 if (patch->partition_of_unity && patch->multiplicative) { 1670 ierr = VecPointwiseMult(patch->patchY[i], patch->patchY[i], patch->patch_dof_weights[i]);CHKERRQ(ierr); 1671 } 1672 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchY[i], patch->localY, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1673 1674 /* Get the matrix on the patch but with only global bcs applied. 1675 * This matrix is then multiplied with the result from the previous solve 1676 * to obtain by how much the residual changes. */ 1677 if (patch->multiplicative) { 1678 ierr = MatMult(multMat, patch->patchY[i], patch->patchX[i]);CHKERRQ(ierr); 1679 ierr = VecScale(patch->patchX[i], -1.0);CHKERRQ(ierr); 1680 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchX[i], patch->localX, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1681 if (!patch->save_operators) {ierr = MatDestroy(&multMat);CHKERRQ(ierr);} 1682 } 1683 } 1684 } 1685 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 1686 /* XXX: should we do this on the global vector? */ 1687 if (patch->partition_of_unity && !patch->multiplicative) { 1688 ierr = VecPointwiseMult(patch->localY, patch->localY, patch->dof_weights);CHKERRQ(ierr); 1689 } 1690 /* Now patch->localY contains the solution of the patch solves, so we need to combine them all. */ 1691 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 1692 ierr = VecGetArray(y, &globalY);CHKERRQ(ierr); 1693 ierr = VecGetArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1694 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1695 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1696 ierr = VecRestoreArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1697 1698 /* Now we need to send the global BC values through */ 1699 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1700 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 1701 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1702 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 1703 for (bc = 0; bc < numBcs; ++bc) { 1704 const PetscInt idx = bcNodes[bc]; 1705 if (idx < n) globalY[idx] = globalX[idx]; 1706 } 1707 1708 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1709 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1710 ierr = VecRestoreArray(y, &globalY);CHKERRQ(ierr); 1711 1712 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 1713 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 1717 static PetscErrorCode PCReset_PATCH(PC pc) 1718 { 1719 PC_PATCH *patch = (PC_PATCH *) pc->data; 1720 PetscInt i; 1721 PetscErrorCode ierr; 1722 1723 PetscFunctionBegin; 1724 /* TODO: Get rid of all these ifs */ 1725 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 1726 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 1727 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 1728 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 1729 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 1730 ierr = PetscSectionDestroy(&patch->bcCounts);CHKERRQ(ierr); 1731 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 1732 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 1733 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 1734 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 1735 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 1736 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 1737 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 1738 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 1739 1740 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 1741 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 1742 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 1743 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 1744 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 1745 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 1746 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 1747 1748 if (patch->bcs) { 1749 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->bcs[i]);CHKERRQ(ierr);} 1750 ierr = PetscFree(patch->bcs);CHKERRQ(ierr); 1751 } 1752 if (patch->multBcs) { 1753 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->multBcs[i]);CHKERRQ(ierr);} 1754 ierr = PetscFree(patch->multBcs);CHKERRQ(ierr); 1755 } 1756 if (patch->ksp) { 1757 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset(patch->ksp[i]);CHKERRQ(ierr);} 1758 } 1759 1760 ierr = VecDestroy(&patch->localX);CHKERRQ(ierr); 1761 ierr = VecDestroy(&patch->localY);CHKERRQ(ierr); 1762 if (patch->patchX) { 1763 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchX[i]);CHKERRQ(ierr);} 1764 ierr = PetscFree(patch->patchX);CHKERRQ(ierr); 1765 } 1766 if (patch->patchY) { 1767 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchY[i]);CHKERRQ(ierr);} 1768 ierr = PetscFree(patch->patchY);CHKERRQ(ierr); 1769 } 1770 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 1771 if (patch->patch_dof_weights) { 1772 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 1773 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 1774 } 1775 if (patch->mat) { 1776 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 1777 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 1778 } 1779 if (patch->multMat) { 1780 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->multMat[i]);CHKERRQ(ierr);} 1781 ierr = PetscFree(patch->multMat);CHKERRQ(ierr); 1782 } 1783 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 1784 if (patch->userIS) { 1785 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 1786 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 1787 } 1788 patch->bs = 0; 1789 patch->cellNodeMap = NULL; 1790 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 1791 1792 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 1793 PetscFunctionReturn(0); 1794 } 1795 1796 static PetscErrorCode PCDestroy_PATCH(PC pc) 1797 { 1798 PC_PATCH *patch = (PC_PATCH *) pc->data; 1799 PetscInt i; 1800 PetscErrorCode ierr; 1801 1802 PetscFunctionBegin; 1803 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 1804 if (patch->ksp) { 1805 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy(&patch->ksp[i]);CHKERRQ(ierr);} 1806 ierr = PetscFree(patch->ksp);CHKERRQ(ierr); 1807 } 1808 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1809 PetscFunctionReturn(0); 1810 } 1811 1812 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 1813 { 1814 PC_PATCH *patch = (PC_PATCH *) pc->data; 1815 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 1816 char sub_mat_type[PETSC_MAX_PATH_LEN]; 1817 const char *prefix; 1818 PetscBool flg, dimflg, codimflg; 1819 MPI_Comm comm; 1820 PetscErrorCode ierr; 1821 1822 PetscFunctionBegin; 1823 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 1824 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 1825 ierr = PetscOptionsHead(PetscOptionsObject, "Vertex-patch Additive Schwarz options");CHKERRQ(ierr); 1826 ierr = PetscOptionsBool("-pc_patch_save_operators", "Store all patch operators for lifetime of PC?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 1827 ierr = PetscOptionsBool("-pc_patch_partition_of_unity", "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 1828 ierr = PetscOptionsBool("-pc_patch_multiplicative", "Gauss-Seidel instead of Jacobi?", "PCPatchSetMultiplicative", patch->multiplicative, &patch->multiplicative, &flg);CHKERRQ(ierr); 1829 ierr = PetscOptionsInt("-pc_patch_construct_dim", "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);CHKERRQ(ierr); 1830 ierr = PetscOptionsInt("-pc_patch_construct_codim", "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);CHKERRQ(ierr); 1831 if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr); 1832 ierr = PetscOptionsEnum("-pc_patch_construct_type", "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 1833 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 1834 ierr = PetscOptionsInt("-pc_patch_vanka_dim", "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 1835 ierr = PetscOptionsInt("-pc_patch_ignore_dim", "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 1836 ierr = PetscOptionsFList("-pc_patch_sub_mat_type", "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);CHKERRQ(ierr); 1837 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 1838 ierr = PetscOptionsBool("-pc_patch_symmetrise_sweep", "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 1839 ierr = PetscOptionsInt("-pc_patch_exclude_subspace", "What subspace (if any) to exclude in construction?", "PCPATCH", patch->exclude_subspace, &patch->exclude_subspace, &flg);CHKERRQ(ierr); 1840 1841 ierr = PetscOptionsBool("-pc_patch_patches_view", "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 1842 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_cells_view", &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 1843 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_points_view", &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 1844 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_section_view", &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 1845 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_sub_mat_view", &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 1846 ierr = PetscOptionsTail();CHKERRQ(ierr); 1847 patch->optionsSet = PETSC_TRUE; 1848 PetscFunctionReturn(0); 1849 } 1850 1851 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 1852 { 1853 PC_PATCH *patch = (PC_PATCH*) pc->data; 1854 KSPConvergedReason reason; 1855 PetscInt i; 1856 PetscErrorCode ierr; 1857 1858 PetscFunctionBegin; 1859 for (i = 0; i < patch->npatch; ++i) { 1860 ierr = KSPSetUp(patch->ksp[i]);CHKERRQ(ierr); 1861 ierr = KSPGetConvergedReason(patch->ksp[i], &reason);CHKERRQ(ierr); 1862 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 1863 } 1864 PetscFunctionReturn(0); 1865 } 1866 1867 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 1868 { 1869 PC_PATCH *patch = (PC_PATCH *) pc->data; 1870 PetscViewer sviewer; 1871 PetscBool isascii; 1872 PetscMPIInt rank; 1873 PetscErrorCode ierr; 1874 1875 PetscFunctionBegin; 1876 /* TODO Redo tabbing with set tbas in new style */ 1877 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 1878 if (!isascii) PetscFunctionReturn(0); 1879 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 1880 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1881 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 1882 if (patch->multiplicative) {ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");CHKERRQ(ierr);} 1883 else {ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr);} 1884 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 1885 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 1886 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 1887 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 1888 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 1889 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 1890 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 1891 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 1892 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 1893 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 1894 ierr = PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");CHKERRQ(ierr); 1895 if (patch->ksp) { 1896 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1897 if (!rank) { 1898 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 1899 ierr = KSPView(patch->ksp[0], sviewer);CHKERRQ(ierr); 1900 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 1901 } 1902 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1903 } else { 1904 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1905 ierr = PetscViewerASCIIPrintf(viewer, "KSP not yet set.\n");CHKERRQ(ierr); 1906 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1907 } 1908 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1909 PetscFunctionReturn(0); 1910 } 1911 1912 /*MC 1913 PCPATCH = "patch" - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 1914 small block additive and multiplicative preconditioners. Block definition is based on topology from 1915 a DM and equation numbering from a PetscSection. 1916 1917 Options Database Keys: 1918 + -pc_patch_cells_view - Views the process local cell numbers for each patch 1919 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 1920 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 1921 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 1922 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 1923 1924 Level: intermediate 1925 1926 .seealso: PCType, PCCreate(), PCSetType() 1927 M*/ 1928 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 1929 { 1930 PC_PATCH *patch; 1931 PetscErrorCode ierr; 1932 1933 PetscFunctionBegin; 1934 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 1935 1936 /* Set some defaults */ 1937 patch->combined = PETSC_FALSE; 1938 patch->save_operators = PETSC_TRUE; 1939 patch->partition_of_unity = PETSC_FALSE; 1940 patch->multiplicative = PETSC_FALSE; 1941 patch->codim = -1; 1942 patch->dim = -1; 1943 patch->exclude_subspace = -1; 1944 patch->vankadim = -1; 1945 patch->ignoredim = -1; 1946 patch->patchconstructop = PCPatchConstruct_Star; 1947 patch->symmetrise_sweep = PETSC_FALSE; 1948 patch->npatch = 0; 1949 patch->userIS = NULL; 1950 patch->optionsSet = PETSC_FALSE; 1951 patch->iterationSet = NULL; 1952 patch->user_patches = PETSC_FALSE; 1953 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 1954 patch->viewPatches = PETSC_FALSE; 1955 patch->viewCells = PETSC_FALSE; 1956 patch->viewPoints = PETSC_FALSE; 1957 patch->viewSection = PETSC_FALSE; 1958 patch->viewMatrix = PETSC_FALSE; 1959 1960 pc->data = (void *) patch; 1961 pc->ops->apply = PCApply_PATCH; 1962 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 1963 pc->ops->setup = PCSetUp_PATCH; 1964 pc->ops->reset = PCReset_PATCH; 1965 pc->ops->destroy = PCDestroy_PATCH; 1966 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 1967 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 1968 pc->ops->view = PCView_PATCH; 1969 pc->ops->applyrichardson = 0; 1970 1971 PetscFunctionReturn(0); 1972 } 1973