xref: /petsc/src/dm/interface/dmgeommodel.c (revision f3d3cd90648576fafae91a24e2611daaef7bcd2e)
1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2 
3 PetscFunctionList DMGeomModelList              = NULL;
4 PetscBool         DMGeomModelRegisterAllCalled = PETSC_FALSE;
5 
6 #if defined(PETSC_HAVE_EGADS)
7 PETSC_EXTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
8 PETSC_EXTERN PetscErrorCode DMSnapToGeomModel_EGADSLite(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
9 #endif
10 
11 static PetscErrorCode DMSnapToGeomModelBall(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
12 {
13   PetscInt val;
14 
15   PetscFunctionBeginUser;
16   PetscCall(DMGetLabelValue(dm, "marker", p, &val));
17   if (val >= 0) {
18     PetscReal norm = 0.;
19 
20     for (PetscInt d = 0; d < dE; ++d) norm += PetscSqr(PetscRealPart(mcoords[d]));
21     norm = PetscSqrtReal(norm);
22     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d] / norm;
23   } else {
24     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
25   }
26   PetscFunctionReturn(PETSC_SUCCESS);
27 }
28 
29 static PetscErrorCode DMSnapToGeomModelCylinder(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
30 {
31   PetscReal gmin[3], gmax[3];
32   PetscInt  val;
33 
34   PetscFunctionBeginUser;
35   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
36   PetscCall(DMGetLabelValue(dm, "generatrix", p, &val));
37   if (val >= 0) {
38     PetscReal norm = 0.;
39 
40     for (PetscInt d = 0; d < dE - 1; ++d) norm += PetscSqr(PetscRealPart(mcoords[d]));
41     norm = PetscSqrtReal(norm);
42     for (PetscInt d = 0; d < dE - 1; ++d) gcoords[d] = mcoords[d] * gmax[0] / norm;
43     gcoords[dE - 1] = mcoords[dE - 1];
44   } else {
45     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
46   }
47   PetscFunctionReturn(PETSC_SUCCESS);
48 }
49 
50 /*@C
51   DMGeomModelRegisterAll - Registers all of the geometry model methods in the `DM` package.
52 
53   Not Collective
54 
55   Level: advanced
56 
57 .seealso: `DM`, `DMGeomModelRegisterDestroy()`
58 @*/
59 PetscErrorCode DMGeomModelRegisterAll(void)
60 {
61   PetscFunctionBegin;
62   if (DMGeomModelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
63   DMGeomModelRegisterAllCalled = PETSC_TRUE;
64   PetscCall(DMGeomModelRegister("ball", DMSnapToGeomModelBall));
65   PetscCall(DMGeomModelRegister("cylinder", DMSnapToGeomModelCylinder));
66 #if defined(PETSC_HAVE_EGADS)
67   // FIXME: Brandon uses DMPlexSnapToGeomModel() here instead
68   PetscCall(DMGeomModelRegister("egads", DMSnapToGeomModel_EGADS));
69 #endif
70   PetscFunctionReturn(PETSC_SUCCESS);
71 }
72 
73 /*@C
74   DMGeomModelRegister -  Adds a geometry model to `DM`
75 
76   Not Collective, No Fortran Support
77 
78   Input Parameters:
79 + sname - name of a new user-defined geometry model
80 - fnc   - geometry model function
81 
82   Example Usage:
83 .vb
84    DMGeomModelRegister("my_geom_model", MySnapToGeomModel);
85 .ve
86 
87   Then, your generator can be chosen with the procedural interface via
88 .vb
89   DMSetGeomModel(dm, "my_geom_model",...)
90 .ve
91   or at runtime via the option
92 .vb
93   -dm_geom_model my_geom_model
94 .ve
95 
96   Level: advanced
97 
98   Note:
99   `DMGeomModelRegister()` may be called multiple times to add several user-defined generators
100 
101 .seealso: `DM`, `DMGeomModelRegisterAll()`, `DMPlexGeomModel()`, `DMGeomModelRegisterDestroy()`
102 @*/
103 PetscErrorCode DMGeomModelRegister(const char sname[], PetscErrorCode (*fnc)(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]))
104 {
105   PetscFunctionBegin;
106   PetscCall(PetscFunctionListAdd(&DMGeomModelList, sname, fnc));
107   PetscFunctionReturn(PETSC_SUCCESS);
108 }
109 
110 extern PetscBool DMGeomModelRegisterAllCalled;
111 
112 PetscErrorCode DMGeomModelRegisterDestroy(void)
113 {
114   PetscFunctionBegin;
115   PetscCall(PetscFunctionListDestroy(&DMGeomModelList));
116   DMGeomModelRegisterAllCalled = PETSC_FALSE;
117   PetscFunctionReturn(PETSC_SUCCESS);
118 }
119 
120 /*@
121   DMSetSnapToGeomModel - Choose a geometry model for this `DM`.
122 
123   Not Collective
124 
125   Input Parameters:
126 + dm   - The `DM` object
127 - name - A geometry model name, or `NULL` for the default
128 
129   Level: intermediate
130 
131 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMSnapToGeomModel()`
132 @*/
133 PetscErrorCode DMSetSnapToGeomModel(DM dm, const char name[])
134 {
135   char      geomname[PETSC_MAX_PATH_LEN];
136   PetscBool flg;
137 
138   PetscFunctionBegin;
139   if (!name && dm->ops->snaptogeommodel) PetscFunctionReturn(PETSC_SUCCESS);
140   PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_geom_model", geomname, sizeof(geomname), &flg));
141   if (flg) name = geomname;
142   if (!name) {
143     PetscObject modelObj;
144 
145     PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", &modelObj));
146     if (modelObj) name = "egads";
147     else {
148       PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", &modelObj));
149       if (modelObj) name = "egads";
150     }
151   }
152   if (!name) PetscFunctionReturn(PETSC_SUCCESS);
153 
154   PetscCall(PetscFunctionListFind(DMGeomModelList, name, &dm->ops->snaptogeommodel));
155   PetscCheck(dm->ops->snaptogeommodel, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Geometry model %s not registered; you may need to add --download-%s to your ./configure options", name, name);
156   PetscFunctionReturn(PETSC_SUCCESS);
157 }
158 
159 /*@
160   DMSnapToGeomModel - Given a coordinate point 'mcoords' on the mesh point 'p', return the closest coordinate point 'gcoords' on the geometry model associated with that point.
161 
162   Not Collective
163 
164   Input Parameters:
165 + dm      - The `DMPLEX` object
166 . p       - The mesh point
167 . dE      - The coordinate dimension
168 - mcoords - A coordinate point lying on the mesh point
169 
170   Output Parameter:
171 . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point
172 
173   Level: intermediate
174 
175   Note:
176   Returns the original coordinates if no geometry model is found.
177 
178   The coordinate dimension may be different from the coordinate dimension of the `dm`, for example if the transformation is extrusion.
179 
180 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMPlexSetRefinementUniform()`
181 @*/
182 PetscErrorCode DMSnapToGeomModel(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
183 {
184   PetscFunctionBegin;
185   if (!dm->ops->snaptogeommodel)
186     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
187   else PetscUseTypeMethod(dm, snaptogeommodel, p, dE, mcoords, gcoords);
188   PetscFunctionReturn(PETSC_SUCCESS);
189 }
190