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:   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: }