Actual source code: plexhpddm.c
1: #include <petsc/private/dmpleximpl.h>
3: /*
4: DMCreateNeumannOverlap - Generates an IS, an unassembled (Neumann) Mat, a setup function, and the corresponding context to be used by PCHPDDM.
6: Input Parameter:
7: . dm - preconditioner context
9: Output Parameters:
10: + ovl - index set of overlapping subdomains
11: . J - unassembled (Neumann) local matrix
12: . setup - function for generating the matrix
13: - setup_ctx - context for setup
15: Options Database Keys:
16: + -dm_plex_view_neumann_original - view the DM without overlap
17: - -dm_plex_view_neumann_overlap - view the DM with overlap as needed by PCHPDDM
19: Level: advanced
21: .seealso: `DMCreate()`, `DM`, `MATIS`, `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()`
22: */
23: PetscErrorCode DMCreateNeumannOverlap_Plex(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx)
24: {
25: DM odm;
26: Mat pJ;
27: PetscSF sf = NULL;
28: PetscSection sec, osec;
29: ISLocalToGlobalMapping l2g;
30: const PetscInt *idxs;
31: PetscInt n, mh;
33: PetscFunctionBegin;
34: *setup = NULL;
35: *setup_ctx = NULL;
36: *ovl = NULL;
37: *J = NULL;
39: /* Overlapped mesh
40: overlap is a little more generous, since it is not computed starting from the owned (Dirichlet) points, but from the locally owned cells */
41: PetscCall(DMPlexDistributeOverlap(dm, 1, &sf, &odm));
42: if (!odm) {
43: PetscCall(PetscSFDestroy(&sf));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /* share discretization */
48: PetscCall(DMGetLocalSection(dm, &sec));
49: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec));
50: PetscCall(PetscSFDistributeSection(sf, sec, NULL, osec));
51: /* what to do here? using both is fine? */
52: PetscCall(DMSetLocalSection(odm, osec));
53: PetscCall(DMCopyDisc(dm, odm));
54: PetscCall(DMPlexGetMaxProjectionHeight(dm, &mh));
55: PetscCall(DMPlexSetMaxProjectionHeight(odm, mh));
56: PetscCall(PetscSectionDestroy(&osec));
58: /* material parameters */
59: {
60: Vec A;
62: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &A));
63: if (A) {
64: DM dmAux, ocdm, odmAux;
65: Vec oA;
67: PetscCall(VecGetDM(A, &dmAux));
68: PetscCall(DMClone(odm, &odmAux));
69: PetscCall(DMGetCoordinateDM(odm, &ocdm));
70: PetscCall(DMSetCoordinateDM(odmAux, ocdm));
71: PetscCall(DMCopyDisc(dmAux, odmAux));
73: PetscCall(DMGetLocalSection(dmAux, &sec));
74: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec));
75: PetscCall(VecCreate(PetscObjectComm((PetscObject)A), &oA));
76: PetscCall(VecSetDM(oA, odmAux));
77: /* TODO: what if these values changes? */
78: PetscCall(DMPlexDistributeField(dmAux, sf, sec, A, osec, oA));
79: PetscCall(DMSetLocalSection(odmAux, osec));
80: PetscCall(PetscSectionDestroy(&osec));
81: PetscCall(DMSetAuxiliaryVec(odm, NULL, 0, 0, oA));
82: PetscCall(VecDestroy(&oA));
83: PetscCall(DMDestroy(&odmAux));
84: }
85: }
86: PetscCall(PetscSFDestroy(&sf));
88: PetscCall(DMViewFromOptions(dm, NULL, "-dm_plex_view_neumann_original"));
89: PetscCall(PetscObjectSetName((PetscObject)odm, "OVL"));
90: PetscCall(DMViewFromOptions(odm, NULL, "-dm_plex_view_neumann_overlap"));
92: /* MATIS for the overlap region
93: the HPDDM interface wants local matrices,
94: we attach the global MATIS to the overlap IS,
95: since we need it to do assembly */
96: PetscCall(DMSetMatType(odm, MATIS));
97: PetscCall(DMCreateMatrix(odm, &pJ));
98: PetscCall(MatISGetLocalMat(pJ, J));
99: PetscCall(PetscObjectReference((PetscObject)*J));
101: /* overlap IS */
102: PetscCall(MatISGetLocalToGlobalMapping(pJ, &l2g, NULL));
103: PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n));
104: PetscCall(ISLocalToGlobalMappingGetIndices(l2g, &idxs));
105: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)odm), n, idxs, PETSC_COPY_VALUES, ovl));
106: PetscCall(ISLocalToGlobalMappingRestoreIndices(l2g, &idxs));
107: PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject)pJ));
108: PetscCall(DMDestroy(&odm));
109: PetscCall(MatDestroy(&pJ));
111: /* special purpose setup function (composed in DMPlexSetSNESLocalFEM) */
112: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup));
113: if (*setup) PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }