Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/libmesh_config.h"
22 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 : #include "libmesh/inf_fe_map.h"
24 : #include "libmesh/inf_fe.h"
25 : #include "libmesh/fe.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/inf_fe_macro.h"
28 : #include "libmesh/libmesh_logging.h"
29 :
30 : #include "libmesh/type_tensor.h"
31 :
32 : namespace libMesh
33 : {
34 :
35 : // ------------------------------------------------------------
36 : // InfFE static class members concerned with coordinate
37 : // mapping
38 :
39 :
40 225020 : Point InfFEMap::map (const unsigned int dim,
41 : const Elem * inf_elem,
42 : const Point & reference_point)
43 : {
44 75359 : libmesh_assert(inf_elem);
45 75359 : libmesh_assert_not_equal_to (dim, 0);
46 :
47 225020 : std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem (inf_elem);
48 :
49 225020 : const Real v (reference_point(dim-1));
50 :
51 : // map in the base face
52 75359 : Point base_point;
53 225020 : switch (dim)
54 : {
55 0 : case 1:
56 0 : base_point = inf_elem->point(0);
57 0 : break;
58 75359 : case 2:
59 : case 3:
60 225020 : base_point = FEMap::map (dim-1, base_elem.get(), reference_point);
61 75359 : break;
62 0 : default:
63 : #ifdef DEBUG
64 0 : libmesh_error_msg("Unknown dim = " << dim);
65 : #endif
66 : break;
67 : }
68 :
69 : // This is the same as the algorithm used below,
70 : // but is more explicit in the actual form in the end.
71 : //
72 : // NOTE: the form used below can be implemented to yield
73 : // more general/flexible mappings, but the current form is
74 : // used e.g. for \p inverse_map() and \p reinit() explicitly.
75 450040 : return (base_point-inf_elem->origin())*2/(1-v)+inf_elem->origin();
76 :
77 : #if 0 // Leave this documented, but w/o unreachable-code warnings
78 : const Order radial_mapping_order (InfFERadial::mapping_order());
79 :
80 : // map in the outer node face not necessary. Simply
81 : // compute the outer_point = base_point + (base_point-origin)
82 : const Point outer_point (base_point * 2. - inf_elem->origin());
83 :
84 : Point p;
85 :
86 : // there are only two mapping shapes in radial direction
87 : p.add_scaled (base_point, eval (v, radial_mapping_order, 0));
88 : p.add_scaled (outer_point, eval (v, radial_mapping_order, 1));
89 :
90 : return p;
91 : #endif
92 74659 : }
93 :
94 :
95 :
96 2415 : Point InfFEMap::inverse_map (const unsigned int dim,
97 : const Elem * inf_elem,
98 : const Point & physical_point,
99 : const Real tolerance,
100 : const bool secure)
101 : {
102 357 : libmesh_assert(inf_elem);
103 357 : libmesh_assert_greater_equal (tolerance, 0.);
104 357 : libmesh_assert(dim > 0);
105 :
106 : // Start logging the map inversion.
107 714 : LOG_SCOPE("inverse_map()", "InfFEMap");
108 :
109 : // The strategy is:
110 : // compute the intersection of the line
111 : // physical_point - origin with the base element,
112 : // find its internal coordinatels using FEMap::inverse_map():
113 : // The radial part can then be computed directly later on.
114 :
115 : // 1.)
116 : // build a base element to do the map inversion in the base face
117 2772 : std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem (inf_elem);
118 :
119 : // The point on the reference element (which we are looking for).
120 : // start with an invalid guess:
121 357 : Point p;
122 2415 : p(dim-1)=-2.;
123 :
124 : // 2.)
125 : // Now find the intersection of a plane represented by the base
126 : // element nodes and the line given by the origin of the infinite
127 : // element and the physical point.
128 357 : Point intersection;
129 :
130 : // the origin of the infinite element
131 2415 : const Point o = inf_elem->origin();
132 :
133 2415 : switch (dim)
134 : {
135 : // unnecessary for 1D
136 0 : case 1:
137 0 : break;
138 :
139 0 : case 2:
140 0 : libmesh_error_msg("ERROR: InfFE::inverse_map is not yet implemented in 2d");
141 :
142 357 : case 3:
143 : {
144 711 : const Point xi ( base_elem->point(1) - base_elem->point(0));
145 357 : const Point eta( base_elem->point(2) - base_elem->point(0));
146 357 : const Point zeta( physical_point - o);
147 :
148 : // normal vector of the base elements plane
149 2061 : Point n(xi.cross(eta));
150 2415 : Real c_factor = (base_elem->point(0) - o)*n/(zeta*n) - 1.;
151 :
152 : // Check whether the above system is ill-posed. It should
153 : // only happen when \p physical_point is not in \p
154 : // inf_elem. In this case, any point that is not in the
155 : // element is a valid answer.
156 2415 : if (libmesh_isinf(c_factor))
157 0 : return p;
158 :
159 : // Compute the intersection with
160 : // {intersection} = {physical_point} + c*({physical_point}-{o}).
161 2415 : intersection.add_scaled(physical_point,1.+c_factor);
162 2415 : intersection.add_scaled(o,-c_factor);
163 :
164 : // For non-planar elements, the intersecting point is obtained via Newton-iteration
165 2415 : if (!base_elem->has_affine_map())
166 : {
167 0 : unsigned int iter_max = 20;
168 :
169 : // the number of shape functions needed for the base_elem
170 1390 : unsigned int n_sf = FE<2,LAGRANGE>::n_shape_functions(base_elem->type(),base_elem->default_order());
171 :
172 : // guess base element coordinates: p=xi,eta,0
173 : // in first iteration, never run with 'secure' to avoid false warnings.
174 1390 : Point ref_point= FEMap::inverse_map(dim-1, base_elem.get(), intersection,
175 0 : tolerance, false);
176 :
177 : // Newton iteration
178 2327 : for(unsigned int it=0; it<iter_max; ++it)
179 : {
180 : // Get the shape function and derivative values at the reference coordinate
181 : // phi.size() == dphi.size()
182 0 : Point dxyz_dxi;
183 0 : Point dxyz_deta;
184 :
185 0 : Point intersection_guess;
186 23270 : for(unsigned int i=0; i<n_sf; ++i)
187 : {
188 :
189 20943 : intersection_guess += base_elem->node_ref(i) * FE<2,LAGRANGE>::shape(base_elem->type(),
190 20943 : base_elem->default_order(),
191 : i,
192 0 : ref_point);
193 :
194 20943 : dxyz_dxi += base_elem->node_ref(i) * FE<2,LAGRANGE>::shape_deriv(base_elem->type(),
195 20943 : base_elem->default_order(),
196 : i,
197 : 0, // d()/dxi
198 0 : ref_point);
199 :
200 20943 : dxyz_deta+= base_elem->node_ref(i) * FE<2,LAGRANGE>::shape_deriv(base_elem->type(),
201 20943 : base_elem->default_order(),
202 : i,
203 : 1, // d()/deta
204 0 : ref_point);
205 : } // for i
206 :
207 0 : TypeVector<Real> F(physical_point + c_factor*(physical_point-o) - intersection_guess);
208 :
209 0 : TypeTensor<Real> J;
210 2327 : J(0,0) = (physical_point-o)(0);
211 2327 : J(0,1) = -dxyz_dxi(0);
212 2327 : J(0,2) = -dxyz_deta(0);
213 2327 : J(1,0) = (physical_point-o)(1);
214 2327 : J(1,1) = -dxyz_dxi(1);
215 2327 : J(1,2) = -dxyz_deta(1);
216 2327 : J(2,0) = (physical_point-o)(2);
217 2327 : J(2,1) = -dxyz_dxi(2);
218 2327 : J(2,2) = -dxyz_deta(2);
219 :
220 : // delta will be the newton step
221 0 : TypeVector<Real> delta;
222 2327 : J.solve(F,delta);
223 :
224 : // check for convergence
225 2327 : Real tol = std::min( TOLERANCE*0.01, TOLERANCE*base_elem->hmax() );
226 2327 : if ( delta.norm() < tol )
227 : {
228 : // newton solver converged, now make sure it converged to a point on the base_elem
229 1390 : if (base_elem->contains_point(intersection_guess,TOLERANCE*0.1))
230 : {
231 1266 : intersection = intersection_guess;
232 : }
233 : //else
234 : // {
235 : // err<<" Error: inverse_map failed!"<<std::endl;
236 : // }
237 0 : break; // break out of 'for it'
238 : }
239 : else
240 : {
241 937 : c_factor -= delta(0);
242 937 : ref_point(0) -= delta(1);
243 937 : ref_point(1) -= delta(2);
244 : }
245 :
246 : }
247 :
248 : }
249 2061 : break;
250 : }
251 :
252 0 : default:
253 0 : libmesh_error_msg("Invalid dim = " << dim);
254 : }
255 :
256 : // 3.)
257 : // Now we have the intersection-point (projection of physical point onto base-element).
258 : // Lets compute its internal coordinates (being p(0) and p(1)):
259 1704 : p= FEMap::inverse_map(dim-1, base_elem.get(), intersection,
260 1065 : tolerance, secure, secure);
261 :
262 : // 4.
263 : // Now that we have the local coordinates in the base,
264 : // we compute the radial distance with Newton iteration.
265 :
266 : // distance from the physical point to the ifem origin
267 2061 : const Real fp_o_dist = Point(o-physical_point).norm();
268 :
269 : // the distance from the intersection on the
270 : // base to the origin
271 2061 : const Real a_dist = Point(o-intersection).norm();
272 :
273 : // element coordinate in radial direction
274 :
275 : // fp_o_dist is at infinity.
276 2415 : if (libmesh_isinf(fp_o_dist))
277 : {
278 0 : p(dim-1)=1;
279 0 : return p;
280 : }
281 :
282 : // when we are somewhere in this element:
283 357 : Real v = 0;
284 :
285 : // For now we're sticking with T_map == CARTESIAN
286 : // if (T_map == CARTESIAN)
287 2415 : v = 1.-2.*a_dist/fp_o_dist;
288 : // else
289 : // libmesh_not_implemented();
290 :
291 :
292 2415 : p(dim-1)=v;
293 : #ifdef DEBUG
294 : // first check whether we are in the reference-element:
295 357 : if ((-1.-1.e-5 < v && v < 1.) || secure)
296 : {
297 357 : const Point check = map (dim, inf_elem, p);
298 357 : const Point diff = physical_point - check;
299 :
300 357 : if (diff.norm() > tolerance)
301 : libmesh_warning("WARNING: diff is " << diff.norm());
302 : }
303 : #endif
304 :
305 2415 : return p;
306 1704 : }
307 :
308 :
309 :
310 0 : void InfFEMap::inverse_map (const unsigned int dim,
311 : const Elem * elem,
312 : const std::vector<Point> & physical_points,
313 : std::vector<Point> & reference_points,
314 : const Real tolerance,
315 : const bool secure)
316 : {
317 : // The number of points to find the
318 : // inverse map of
319 0 : const std::size_t n_points = physical_points.size();
320 :
321 : // Resize the vector to hold the points
322 : // on the reference element
323 0 : reference_points.resize(n_points);
324 :
325 : // Find the coordinates on the reference
326 : // element of each point in physical space
327 0 : for (unsigned int p=0; p<n_points; p++)
328 0 : reference_points[p] =
329 0 : inverse_map (dim, elem, physical_points[p], tolerance, secure);
330 0 : }
331 :
332 :
333 : } // namespace libMesh
334 :
335 :
336 : #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
