xref: /petsc/src/dm/impls/plex/transform/impls/filter/plextrfilter.c (revision 31d78bcd2b98084dc1368b20eb1129c8b9fb39fe)
1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2 
DMPlexTransformView_Filter(DMPlexTransform tr,PetscViewer viewer)3 static PetscErrorCode DMPlexTransformView_Filter(DMPlexTransform tr, PetscViewer viewer)
4 {
5   PetscBool isascii;
6 
7   PetscFunctionBegin;
8   PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
9   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
10   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
11   if (isascii) {
12     const char *name;
13 
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 }
21 
DMPlexTransformSetUp_Filter(DMPlexTransform tr)22 static PetscErrorCode DMPlexTransformSetUp_Filter(DMPlexTransform tr)
23 {
24   DM       dm;
25   DMLabel  active;
26   PetscInt Nc;
27 
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;
35 
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;
45 
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 }
60 
DMPlexTransformDestroy_Filter(DMPlexTransform tr)61 static PetscErrorCode DMPlexTransformDestroy_Filter(DMPlexTransform tr)
62 {
63   DMPlexTransform_Filter *f = (DMPlexTransform_Filter *)tr->data;
64 
65   PetscFunctionBegin;
66   PetscCall(PetscFree(f));
67   PetscFunctionReturn(PETSC_SUCCESS);
68 }
69 
DMPlexTransformCellTransform_Filter(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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;
75 
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 }
91 
DMPlexTransformInitialize_Filter(DMPlexTransform tr)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 }
104 
DMPlexTransformCreate_Filter(DMPlexTransform tr)105 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform tr)
106 {
107   DMPlexTransform_Filter *f;
108 
109   PetscFunctionBegin;
110   PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
111   PetscCall(PetscNew(&f));
112   tr->data = f;
113 
114   PetscCall(DMPlexTransformInitialize_Filter(tr));
115   PetscFunctionReturn(PETSC_SUCCESS);
116 }
117