1: #include <petsc/private/dmplextransformimpl.h>

  3: static PetscErrorCode DMPlexTransformView_Filter(DMPlexTransform tr, PetscViewer viewer)
  4: {
  5:   PetscBool isascii;

  7:   PetscFunctionBegin;
 10:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
 11:   if (isascii) {
 12:     const char *name;

 14:     PetscCall(PetscObjectGetName((PetscObject)tr, &name));
 15:     PetscCall(PetscViewerASCIIPrintf(viewer, "Filter transformation %s\n", name ? name : ""));
 16:   } else {
 17:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
 18:   }
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

 22: static PetscErrorCode DMPlexTransformSetUp_Filter(DMPlexTransform tr)
 23: {
 24:   DM       dm;
 25:   DMLabel  active;
 26:   PetscInt Nc;

 28:   PetscFunctionBegin;
 29:   PetscCall(DMPlexTransformGetDM(tr, &dm));
 30:   PetscCall(DMPlexTransformGetActive(tr, &active));
 31:   if (active) {
 32:     IS              filterIS;
 33:     const PetscInt *filterCells;
 34:     PetscInt        c;

 36:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Filter Type", &tr->trType));
 37:     PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &filterIS));
 38:     PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc));
 39:     if (filterIS) PetscCall(ISGetIndices(filterIS, &filterCells));
 40:     for (c = 0; c < Nc; ++c) {
 41:       const PetscInt cell    = filterCells[c];
 42:       PetscInt      *closure = NULL;
 43:       DMPolytopeType ct;
 44:       PetscInt       Ncl, cl;

 46:       PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
 47:       for (cl = 0; cl < Ncl * 2; cl += 2) {
 48:         PetscCall(DMPlexGetCellType(dm, closure[cl], &ct));
 49:         PetscCall(DMLabelSetValue(tr->trType, closure[cl], ct));
 50:       }
 51:       PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
 52:     }
 53:     if (filterIS) {
 54:       PetscCall(ISRestoreIndices(filterIS, &filterCells));
 55:       PetscCall(ISDestroy(&filterIS));
 56:     }
 57:   }
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 61: static PetscErrorCode DMPlexTransformDestroy_Filter(DMPlexTransform tr)
 62: {
 63:   DMPlexTransform_Filter *f = (DMPlexTransform_Filter *)tr->data;

 65:   PetscFunctionBegin;
 66:   PetscCall(PetscFree(f));
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: static PetscErrorCode DMPlexTransformCellTransform_Filter(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
 71: {
 72:   PetscFunctionBeginHot;
 73:   if (tr->trType && p >= 0) {
 74:     PetscInt val;

 76:     PetscCall(DMLabelGetValue(tr->trType, p, &val));
 77:     if (val >= 0) {
 78:       if (rt) *rt = val;
 79:       PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
 80:       PetscFunctionReturn(PETSC_SUCCESS);
 81:     }
 82:   }
 83:   if (rt) *rt = -1;
 84:   *Nt     = 0;
 85:   *target = NULL;
 86:   *size   = NULL;
 87:   *cone   = NULL;
 88:   *ornt   = NULL;
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode DMPlexTransformInitialize_Filter(DMPlexTransform tr)
 93: {
 94:   PetscFunctionBegin;
 95:   tr->ops->view                  = DMPlexTransformView_Filter;
 96:   tr->ops->setup                 = DMPlexTransformSetUp_Filter;
 97:   tr->ops->destroy               = DMPlexTransformDestroy_Filter;
 98:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Internal;
 99:   tr->ops->celltransform         = DMPlexTransformCellTransform_Filter;
100:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientationIdentity;
101:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform tr)
106: {
107:   DMPlexTransform_Filter *f;

109:   PetscFunctionBegin;
111:   PetscCall(PetscNew(&f));
112:   tr->data = f;

114:   PetscCall(DMPlexTransformInitialize_Filter(tr));
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }