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:     DM  dmAux;
 62:     Vec A;

 64:     PetscObjectQuery((PetscObject)dm, "dmAux", (PetscObject*)&dmAux);
 65:     PetscObjectQuery((PetscObject)dm, "A", (PetscObject*)&A);
 66:     if (dmAux) {
 67:       DM ocdm, odmAux;

 69:       DMClone(odm, &odmAux);
 70:       DMGetCoordinateDM(odm, &ocdm);
 71:       DMSetCoordinateDM(odmAux, ocdm);
 72:       DMCopyDisc(dmAux, odmAux);
 73:       if (A) {
 74:         Vec oA;

 76:         DMGetLocalSection(dmAux, &sec);
 77:         PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec);
 78:         VecCreate(PetscObjectComm((PetscObject)A), &oA);
 79:         /* TODO: what if these values changes? */
 80:         DMPlexDistributeField(dmAux, sf, sec, A, osec, oA);
 81:         DMSetLocalSection(odmAux, osec);
 82:         PetscSectionDestroy(&osec);
 83:         PetscObjectCompose((PetscObject)odm, "A", (PetscObject)oA);
 84:         VecDestroy(&oA);
 85:       }
 86:       PetscObjectCompose((PetscObject)odm, "dmAux", (PetscObject)odmAux);
 87:       DMDestroy(&odmAux);
 88:     }
 89:   }
 90:   PetscSFDestroy(&sf);

 92:   DMViewFromOptions(dm, NULL, "-dm_plex_view_neumann_original");
 93:   PetscObjectSetName((PetscObject)odm, "OVL");
 94:   DMViewFromOptions(odm, NULL, "-dm_plex_view_neumann_overlap");

 96:   /* MATIS for the overlap region
 97:      the HPDDM interface wants local matrices,
 98:      we attach the global MATIS to the overlap IS,
 99:      since we need it to do assembly */
100:   DMSetMatType(odm, MATIS);
101:   DMCreateMatrix(odm, &pJ);
102:   MatISGetLocalMat(pJ, J);
103:   PetscObjectReference((PetscObject)*J);

105:   /* overlap IS */
106:   MatGetLocalToGlobalMapping(pJ, &l2g, NULL);
107:   ISLocalToGlobalMappingGetSize(l2g, &n);
108:   ISLocalToGlobalMappingGetIndices(l2g, &idxs);
109:   ISCreateGeneral(PetscObjectComm((PetscObject)odm), n, idxs, PETSC_COPY_VALUES, ovl);
110:   ISLocalToGlobalMappingRestoreIndices(l2g, &idxs);
111:   PetscObjectCompose((PetscObject)*ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject)pJ);
112:   DMDestroy(&odm);
113:   MatDestroy(&pJ);

115:   /* special purpose setup function (composed in DMPlexSetSNESLocalFEM) */
116:   PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup);
117:   if (*setup) {
118:     PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm);
119:   }
120:   return(0);
121: }