Actual source code: plextrfilter.c

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