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;

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

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

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

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

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

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

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

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

108:   PetscFunctionBegin;
110:   PetscCall(PetscNew(&f));
111:   tr->redFactor = 1.0;
112:   tr->data      = f;

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