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