1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2
DMPlexTransformSetUp_1D(DMPlexTransform tr)3 static PetscErrorCode DMPlexTransformSetUp_1D(DMPlexTransform tr)
4 {
5 DM dm;
6 DMLabel active;
7 PetscInt pStart, pEnd, p;
8
9 PetscFunctionBegin;
10 PetscCall(DMPlexTransformGetDM(tr, &dm));
11 PetscCall(DMPlexTransformGetActive(tr, &active));
12 PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use 1D algorithm");
13 /* Calculate refineType for each cell */
14 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
15 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
16 for (p = pStart; p < pEnd; ++p) {
17 DMLabel trType = tr->trType;
18 DMPolytopeType ct;
19 PetscInt val;
20
21 PetscCall(DMPlexGetCellType(dm, p, &ct));
22 switch (ct) {
23 case DM_POLYTOPE_POINT:
24 PetscCall(DMLabelSetValue(trType, p, 0));
25 break;
26 case DM_POLYTOPE_SEGMENT:
27 case DM_POLYTOPE_POINT_PRISM_TENSOR:
28 PetscCall(DMLabelGetValue(active, p, &val));
29 if (val == 1) PetscCall(DMLabelSetValue(trType, p, val));
30 else PetscCall(DMLabelSetValue(trType, p, 2));
31 break;
32 default:
33 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
34 }
35 }
36 PetscFunctionReturn(PETSC_SUCCESS);
37 }
38
DMPlexTransformGetSubcellOrientation_1D(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)39 static PetscErrorCode DMPlexTransformGetSubcellOrientation_1D(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
40 {
41 PetscInt rt;
42
43 PetscFunctionBeginHot;
44 PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
45 *rnew = r;
46 *onew = o;
47 switch (rt) {
48 case 1:
49 PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
50 break;
51 default:
52 PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
53 }
54 PetscFunctionReturn(PETSC_SUCCESS);
55 }
56
DMPlexTransformCellTransform_1D(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])57 static PetscErrorCode DMPlexTransformCellTransform_1D(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
58 {
59 DMLabel trType = tr->trType;
60 PetscInt val;
61
62 PetscFunctionBeginHot;
63 PetscCheck(p >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
64 PetscCall(DMLabelGetValue(trType, p, &val));
65 if (rt) *rt = val;
66 switch (source) {
67 case DM_POLYTOPE_POINT:
68 PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
69 break;
70 case DM_POLYTOPE_POINT_PRISM_TENSOR:
71 case DM_POLYTOPE_SEGMENT:
72 if (val == 1) PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
73 else PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
74 break;
75 default:
76 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
77 }
78 PetscFunctionReturn(PETSC_SUCCESS);
79 }
80
DMPlexTransformSetFromOptions_1D(DMPlexTransform tr,PetscOptionItems PetscOptionsObject)81 static PetscErrorCode DMPlexTransformSetFromOptions_1D(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
82 {
83 PetscInt cells[256], n = 256, i;
84 PetscBool flg;
85
86 PetscFunctionBegin;
87 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
88 PetscCall(PetscOptionsIntArray("-dm_plex_transform_1d_ref_cell", "Mark cells for refinement", "", cells, &n, &flg));
89 if (flg) {
90 DMLabel active;
91
92 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active));
93 for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE));
94 PetscCall(DMPlexTransformSetActive(tr, active));
95 PetscCall(DMLabelDestroy(&active));
96 }
97 PetscOptionsHeadEnd();
98 PetscFunctionReturn(PETSC_SUCCESS);
99 }
100
DMPlexTransformView_1D(DMPlexTransform tr,PetscViewer viewer)101 static PetscErrorCode DMPlexTransformView_1D(DMPlexTransform tr, PetscViewer viewer)
102 {
103 PetscBool isascii;
104
105 PetscFunctionBegin;
106 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
107 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
108 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
109 if (isascii) {
110 PetscViewerFormat format;
111 const char *name;
112
113 PetscCall(PetscObjectGetName((PetscObject)tr, &name));
114 PetscCall(PetscViewerASCIIPrintf(viewer, "1D refinement %s\n", name ? name : ""));
115 PetscCall(PetscViewerGetFormat(viewer, &format));
116 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer));
117 } else {
118 SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
119 }
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
DMPlexTransformDestroy_1D(DMPlexTransform tr)123 static PetscErrorCode DMPlexTransformDestroy_1D(DMPlexTransform tr)
124 {
125 PetscFunctionBegin;
126 PetscCall(PetscFree(tr->data));
127 PetscFunctionReturn(PETSC_SUCCESS);
128 }
129
DMPlexTransformInitialize_1D(DMPlexTransform tr)130 static PetscErrorCode DMPlexTransformInitialize_1D(DMPlexTransform tr)
131 {
132 PetscFunctionBegin;
133 tr->ops->view = DMPlexTransformView_1D;
134 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_1D;
135 tr->ops->setup = DMPlexTransformSetUp_1D;
136 tr->ops->destroy = DMPlexTransformDestroy_1D;
137 tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
138 tr->ops->celltransform = DMPlexTransformCellTransform_1D;
139 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_1D;
140 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
DMPlexTransformCreate_1D(DMPlexTransform tr)144 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform tr)
145 {
146 DMPlexRefine_1D *f;
147
148 PetscFunctionBegin;
149 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
150 PetscCall(PetscNew(&f));
151 tr->data = f;
152
153 PetscCall(DMPlexTransformInitialize_1D(tr));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156