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 : // C++ includes
21 :
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_config.h"
25 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
26 : #include "libmesh/inf_fe.h"
27 : #include "libmesh/inf_fe_macro.h"
28 : #include "libmesh/quadrature.h"
29 : #include "libmesh/enum_quadrature_type.h"
30 : #include "libmesh/elem.h"
31 :
32 : namespace libMesh
33 : {
34 :
35 : // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
36 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
37 5 : void InfFE<Dim,T_radial,T_base>::reinit(const Elem * inf_elem,
38 : const unsigned int s,
39 : const Real /*tolerance*/,
40 : const std::vector<Point> * const pts,
41 : const std::vector<Real> * const weights)
42 : {
43 5 : if (weights != nullptr)
44 0 : libmesh_not_implemented_msg("ERROR: User-specified weights for infinite elements are not implemented!");
45 :
46 5 : if (pts != nullptr)
47 0 : libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements are not implemented!");
48 :
49 : // We don't do this for 1D elements!
50 0 : libmesh_assert_not_equal_to (Dim, 1);
51 :
52 2 : libmesh_assert(inf_elem);
53 2 : libmesh_assert(qrule);
54 :
55 : // Build the side of interest
56 5 : const std::unique_ptr<const Elem> side(inf_elem->build_side_ptr(s));
57 :
58 : // set the element and type
59 5 : this->_elem = inf_elem;
60 5 : this->_elem_type = inf_elem->type();
61 :
62 : // eventually initialize radial quadrature rule
63 2 : bool radial_qrule_initialized = false;
64 :
65 : // if we are working on the base-side, the radial function is constant.
66 : // With this, we ensure that at least for base elements we reinitialize all quantities
67 : // when we enter for the first time.
68 5 : if (s == 0)
69 5 : current_fe_type.radial_order = 0;
70 : else
71 : /**
72 : * After the recent larger changes, this case was not tested.
73 : * It might work, but maybe it gives wrong results.
74 : */
75 0 : libmesh_not_implemented();
76 :
77 5 : if (current_fe_type.radial_order != fe_type.radial_order)
78 : {
79 2 : if (s > 0)
80 : {
81 0 : current_fe_type.radial_order = fe_type.radial_order;
82 0 : radial_qrule->init(EDGE2, inf_elem->p_level());
83 : }
84 : else
85 : {
86 : // build a new 0-dimensional quadrature-rule:
87 8 : radial_qrule=QBase::build(QGAUSS, 0, fe_type.radial_order);
88 5 : radial_qrule->init(NODEELEM, 0, /*simple_type_only=*/true);
89 :
90 : //the base_qrule is set up with dim-1, but apparently we need dim, so we replace it:
91 8 : base_qrule=QBase::build(qrule->type(), side->dim(), qrule->get_order());
92 :
93 5 : unsigned int side_p_level = inf_elem->p_level();
94 7 : if (inf_elem->neighbor_ptr(s) != nullptr)
95 7 : side_p_level = std::max(side_p_level, inf_elem->neighbor_ptr(s)->p_level());
96 5 : base_qrule->init(*side, side_p_level);
97 : }
98 2 : radial_qrule_initialized = true;
99 : }
100 :
101 : // Initialize the face shape functions
102 8 : if (this->get_type() != inf_elem->type() ||
103 8 : base_fe->shapes_need_reinit() ||
104 : radial_qrule_initialized)
105 7 : this->init_face_shape_functions (qrule->get_points(), side.get());
106 :
107 : // The reinit() function computes all what we want except for
108 : // - normal, tangents: They are not considered
109 : // This is done below:
110 5 : compute_face_functions();
111 5 : }
112 :
113 :
114 :
115 : // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
116 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
117 0 : void InfFE<Dim,T_radial,T_base>::edge_reinit(const Elem *,
118 : const unsigned int,
119 : const Real,
120 : const std::vector<Point> * const pts,
121 : const std::vector<Real> * const /*weights*/)
122 : {
123 : // We don't do this for 1D elements!
124 : //libmesh_assert_not_equal_to (Dim, 1);
125 0 : libmesh_not_implemented_msg("ERROR: Edge conditions for infinite elements not implemented!");
126 :
127 : if (pts != nullptr)
128 : libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements not implemented!");
129 : }
130 :
131 :
132 :
133 :
134 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
135 5 : void InfFE<Dim,T_radial,T_base>::init_face_shape_functions(const std::vector<Point> &,
136 : const Elem * inf_side)
137 : {
138 2 : libmesh_assert(inf_side);
139 :
140 : // Currently, this makes only sense in 3-D!
141 0 : libmesh_assert_equal_to (Dim, 3);
142 :
143 : // Initialize the radial shape functions (in particular som)
144 5 : this->init_radial_shape_functions(inf_side);
145 :
146 : // Initialize the base shape functions
147 5 : if (inf_side->infinite())
148 0 : this->update_base_elem(inf_side);
149 : else
150 : // in this case, I need the 2D base
151 5 : this->update_base_elem(inf_side->interior_parent());
152 :
153 : // Initialize the base quadrature rule
154 7 : base_qrule->init(*base_elem, inf_side->p_level());
155 :
156 : // base_fe still corresponds to the (dim-1)-dimensional base of the InfFE object,
157 : // so update the fe_base.
158 5 : if (inf_side->infinite())
159 : {
160 0 : base_fe = FEBase::build(Dim-2, this->fe_type);
161 0 : base_fe->attach_quadrature_rule(base_qrule.get());
162 : }
163 : else
164 : {
165 8 : base_fe = FEBase::build(Dim-1, this->fe_type);
166 7 : base_fe->attach_quadrature_rule(base_qrule.get());
167 : }
168 :
169 5 : if (this->calculate_map || this->calculate_map_scaled)
170 : {
171 : //before initializing, we should say what to compute:
172 2 : base_fe->_fe_map->get_xyz();
173 2 : base_fe->_fe_map->get_JxW();
174 : }
175 :
176 5 : if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
177 : // initialize the shape functions on the base
178 9 : base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
179 : base_elem.get());
180 :
181 : // the number of quadrature points
182 2 : const unsigned int n_radial_qp = radial_qrule->n_points();
183 2 : const unsigned int n_base_qp = base_qrule->n_points();
184 5 : const unsigned int n_total_qp = n_radial_qp * n_base_qp;
185 :
186 : #ifdef DEBUG
187 2 : if (som.size() > 0)
188 2 : libmesh_assert_equal_to(n_radial_qp, som.size());
189 : // when evaluating the base side, there should be only one radial point.
190 2 : if (!inf_side->infinite())
191 2 : libmesh_assert_equal_to (n_radial_qp, 1);
192 : #endif
193 :
194 : // the quadrature weights
195 5 : _total_qrule_weights.resize(n_total_qp);
196 9 : std::vector<Point> qp(n_total_qp);
197 :
198 : // quadrature rule weights
199 : if (Dim < 3)
200 : {
201 : // the quadrature points must be assembled differently for lower dims.
202 0 : libmesh_not_implemented();
203 : }
204 : else
205 : {
206 2 : const std::vector<Real> & radial_qw = radial_qrule->get_weights();
207 2 : const std::vector<Real> & base_qw = base_qrule->get_weights();
208 2 : const std::vector<Point> & radial_qp = radial_qrule->get_points();
209 2 : const std::vector<Point> & base_qp = base_qrule->get_points();
210 :
211 2 : libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
212 2 : libmesh_assert_equal_to (base_qw.size(), n_base_qp);
213 :
214 10 : for (unsigned int rp=0; rp<n_radial_qp; rp++)
215 40 : for (unsigned int bp=0; bp<n_base_qp; bp++)
216 : {
217 63 : _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
218 : // initialize the quadrature-points for the 2D side element
219 : // - either the base element or it has a 1D base + radial direction.
220 35 : if (inf_side->infinite())
221 0 : qp[bp + rp*n_base_qp]=Point(base_qp[bp](0),
222 : 0.,
223 0 : radial_qp[rp](0));
224 : else
225 49 : qp[bp + rp*n_base_qp]=Point(base_qp[bp](0),
226 28 : base_qp[bp](1),
227 : -1.);
228 : }
229 : }
230 :
231 5 : this->reinit(inf_side->interior_parent(), &qp);
232 :
233 5 : }
234 :
235 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
236 5 : void InfFE<Dim,T_radial,T_base>::compute_face_functions()
237 : {
238 :
239 5 : if (!calculate_map && !calculate_map_scaled)
240 0 : return; // we didn't ask for any quantity computed here.
241 :
242 4 : const unsigned int n_qp = cast_int<unsigned int>(_total_qrule_weights.size());
243 5 : this->normals.resize(n_qp);
244 :
245 : if (Dim > 1)
246 : {
247 5 : this->tangents.resize(n_qp);
248 40 : for (unsigned int p=0; p<n_qp; ++p)
249 49 : this->tangents[p].resize(LIBMESH_DIM-1);
250 : }
251 : else
252 : {
253 0 : libMesh::err << "tangents have no sense in 1-dimensional elements!"<<std::endl;
254 0 : libmesh_error_msg("Exiting...");
255 : }
256 :
257 : // the dimension of base indicates which side we have:
258 : // if base_dim == Dim -1 : base
259 : // base_dim == Dim -2 : one of the other sides.
260 5 : unsigned int base_dim =base_fe->dim;
261 : // If we have no quadrature points, there's nothing else to do
262 5 : if (!n_qp)
263 0 : return;
264 :
265 : switch(Dim)
266 : {
267 0 : case 1:
268 : case 2:
269 : {
270 0 : libmesh_not_implemented();
271 : break;
272 : }
273 3 : case 3:
274 : {
275 : // Below, we assume a 2D base, i.e. we compute the side s=0.
276 5 : if (base_dim==Dim-1)
277 40 : for (unsigned int p=0; p<n_qp; ++p)
278 : {
279 : //
280 : // seeking dxyzdx, dxyzdeta means to compute
281 : // / dx/dxi dy/dxi dz/dxi \.
282 : // J^-1= | |
283 : // \ dx/deta dy/deta dz/deta /.
284 : // which is the psudo-inverse of J, i.e.
285 : //
286 : // J^-1 = (J^T J)^-1 J^T
287 : //
288 : // where J^T T is the 2x2 matrix 'g' used to compute the
289 : // Jacobian determinant; thus
290 : //
291 : // J^-1 = ________1________ / g22 -g21 \ / dxi/dx dxi/dy dxi/dz \.
292 : // g11*g22 - g21*g12 \-g12 g11 / \ deta/dx deta/dy deta/dz /.
293 35 : const std::vector<Real> & base_dxidx = base_fe->get_dxidx();
294 35 : const std::vector<Real> & base_dxidy = base_fe->get_dxidy();
295 35 : const std::vector<Real> & base_dxidz = base_fe->get_dxidz();
296 35 : const std::vector<Real> & base_detadx = base_fe->get_detadx();
297 35 : const std::vector<Real> & base_detady = base_fe->get_detady();
298 35 : const std::vector<Real> & base_detadz = base_fe->get_detadz();
299 :
300 35 : const Real g11 = (base_dxidx[p]*base_dxidx[p] +
301 35 : base_dxidy[p]*base_dxidy[p] +
302 35 : base_dxidz[p]*base_dxidz[p]);
303 35 : const Real g12 = (base_dxidx[p]*base_detadx[p] +
304 35 : base_dxidy[p]*base_detady[p] +
305 35 : base_dxidz[p]*base_detadz[p]);
306 14 : const Real g21 = g12;
307 63 : const Real g22 = (base_detadx[p]*base_detadx[p] +
308 35 : base_detady[p]*base_detady[p] +
309 35 : base_detadz[p]*base_detadz[p]);
310 :
311 35 : const Real det = (g11*g22 - g12*g21);
312 :
313 49 : Point dxyzdxi_map((g22*base_dxidx[p]-g21*base_detadx[p])/det,
314 35 : (g22*base_dxidy[p]-g21*base_detady[p])/det,
315 35 : (g22*base_dxidz[p]-g21*base_detadz[p])/det);
316 :
317 49 : Point dxyzdeta_map((g11*base_detadx[p] - g12*base_dxidx[p])/det,
318 35 : (g11*base_detady[p] - g12*base_dxidy[p])/det,
319 35 : (g11*base_detadz[p] - g12*base_dxidz[p])/det);
320 :
321 35 : this->tangents[p][0] = dxyzdxi_map.unit();
322 :
323 35 : this->tangents[p][1] = (dxyzdeta_map - (dxyzdeta_map*tangents[p][0])*tangents[p][0] ).unit();
324 :
325 35 : this->normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
326 : // recompute JxW using the 2D Jacobian:
327 : // Since we are at the base, there is no difference between scaled and unscaled jacobian
328 35 : if (calculate_jxw)
329 0 : this->JxW[p] = _total_qrule_weights[p]/std::sqrt(det);
330 :
331 35 : if (calculate_map_scaled)
332 63 : this->JxWxdecay[p] = _total_qrule_weights[p]/std::sqrt(det);
333 :
334 : }
335 0 : else if (base_dim == Dim -2)
336 : {
337 0 : libmesh_not_implemented();
338 : }
339 : else
340 : {
341 : // in this case something went completely wrong.
342 0 : libmesh_not_implemented();
343 : }
344 2 : break;
345 : }
346 : default:
347 : libmesh_error_msg("Unsupported dim = " << dim);
348 : }
349 :
350 : }
351 :
352 :
353 :
354 : // Explicit instantiations - doesn't make sense in 1D, but as
355 : // long as we only return errors, we are fine... ;-)
356 : //#include "libmesh/inf_fe_instantiate_1D.h"
357 : //#include "libmesh/inf_fe_instantiate_2D.h"
358 : //#include "libmesh/inf_fe_instantiate_3D.h"
359 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
360 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
361 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
362 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
363 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
364 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
365 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
366 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
367 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
368 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_face_functions());
369 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_face_functions());
370 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_face_functions());
371 :
372 : } // namespace libMesh
373 :
374 : #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
