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