Actual source code: plexhpddm.c
1: #include <petsc/private/dmpleximpl.h>
3: /*@C
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: *setup = NULL;
34: *setup_ctx = NULL;
35: *ovl = NULL;
36: *J = NULL;
38: /* Overlapped mesh
39: overlap is a little more generous, since it is not computed starting from the owned (Dirichlet) points, but from the locally owned cells */
40: DMPlexDistributeOverlap(dm, 1, &sf, &odm);
41: if (!odm) {
42: PetscSFDestroy(&sf);
43: return 0;
44: }
46: /* share discretization */
47: DMGetLocalSection(dm, &sec);
48: PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec);
49: PetscSFDistributeSection(sf, sec, NULL, osec);
50: /* what to do here? using both is fine? */
51: DMSetLocalSection(odm, osec);
52: DMCopyDisc(dm, odm);
53: DMPlexGetMaxProjectionHeight(dm, &mh);
54: DMPlexSetMaxProjectionHeight(odm, mh);
55: PetscSectionDestroy(&osec);
57: /* material parameters */
58: {
59: Vec A;
61: DMGetAuxiliaryVec(dm, NULL, 0, 0, &A);
62: if (A) {
63: DM dmAux, ocdm, odmAux;
64: Vec oA;
66: VecGetDM(A, &dmAux);
67: DMClone(odm, &odmAux);
68: DMGetCoordinateDM(odm, &ocdm);
69: DMSetCoordinateDM(odmAux, ocdm);
70: DMCopyDisc(dmAux, odmAux);
72: DMGetLocalSection(dmAux, &sec);
73: PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec);
74: VecCreate(PetscObjectComm((PetscObject)A), &oA);
75: VecSetDM(oA, odmAux);
76: /* TODO: what if these values changes? */
77: DMPlexDistributeField(dmAux, sf, sec, A, osec, oA);
78: DMSetLocalSection(odmAux, osec);
79: PetscSectionDestroy(&osec);
80: DMSetAuxiliaryVec(odm, NULL, 0, 0, oA);
81: VecDestroy(&oA);
82: DMDestroy(&odmAux);
83: }
84: }
85: PetscSFDestroy(&sf);
87: DMViewFromOptions(dm, NULL, "-dm_plex_view_neumann_original");
88: PetscObjectSetName((PetscObject)odm, "OVL");
89: DMViewFromOptions(odm, NULL, "-dm_plex_view_neumann_overlap");
91: /* MATIS for the overlap region
92: the HPDDM interface wants local matrices,
93: we attach the global MATIS to the overlap IS,
94: since we need it to do assembly */
95: DMSetMatType(odm, MATIS);
96: DMCreateMatrix(odm, &pJ);
97: MatISGetLocalMat(pJ, J);
98: PetscObjectReference((PetscObject)*J);
100: /* overlap IS */
101: MatISGetLocalToGlobalMapping(pJ, &l2g, NULL);
102: ISLocalToGlobalMappingGetSize(l2g, &n);
103: ISLocalToGlobalMappingGetIndices(l2g, &idxs);
104: ISCreateGeneral(PetscObjectComm((PetscObject)odm), n, idxs, PETSC_COPY_VALUES, ovl);
105: ISLocalToGlobalMappingRestoreIndices(l2g, &idxs);
106: PetscObjectCompose((PetscObject)*ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject)pJ);
107: DMDestroy(&odm);
108: MatDestroy(&pJ);
110: /* special purpose setup function (composed in DMPlexSetSNESLocalFEM) */
111: PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup);
112: if (*setup) PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm);
113: return 0;
114: }