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 : // Local includes
19 : #include "libmesh/libmesh_config.h"
20 :
21 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 :
23 :
24 : // C++ includes
25 : #include <algorithm> // for std::min, std::max
26 :
27 : // Local includes cont'd
28 : #include "libmesh/cell_inf_hex.h"
29 : #include "libmesh/cell_inf_hex8.h"
30 : #include "libmesh/face_quad4.h"
31 : #include "libmesh/face_inf_quad4.h"
32 : #include "libmesh/fe_type.h"
33 : #include "libmesh/fe_interface.h"
34 : #include "libmesh/inf_fe_map.h"
35 : #include "libmesh/enum_elem_quality.h"
36 :
37 : namespace libMesh
38 : {
39 :
40 :
41 :
42 : // ------------------------------------------------------------
43 : // InfHex class static member initializations
44 : const int InfHex::num_sides;
45 : const int InfHex::num_edges;
46 : const int InfHex::num_children;
47 :
48 : const Real InfHex::_master_points[18][3] =
49 : {
50 : {-1, -1, 0},
51 : {1, -1, 0},
52 : {1, 1, 0},
53 : {-1, 1, 0},
54 : {-1, -1, 1},
55 : {1, -1, 1},
56 : {1, 1, 1},
57 : {-1, 1, 1},
58 : {0, -1, 0},
59 : {1, 0, 0},
60 : {0, 1, 0},
61 : {-1, 0, 0},
62 : {0, -1, 1},
63 : {1, 0, 1},
64 : {0, 1, 1},
65 : {-1, 0, 1},
66 : {0, 0, 0},
67 : {0, 0, 1}
68 : };
69 :
70 : const unsigned int InfHex::edge_sides_map[8][2] =
71 : {
72 : {0, 1}, // Edge 0
73 : {0, 2}, // Edge 1
74 : {0, 3}, // Edge 2
75 : {0, 4}, // Edge 3
76 : {1, 4}, // Edge 4
77 : {1, 2}, // Edge 5
78 : {2, 3}, // Edge 6
79 : {3, 4} // Edge 7
80 : };
81 :
82 : const unsigned int InfHex::adjacent_edges_map[/*num_vertices*/8][/*max_adjacent_edges*/3] =
83 : {
84 : {0, 3, 4}, // Edges adjacent to node 0
85 : {0, 1, 5}, // Edges adjacent to node 1
86 : {1, 2, 6}, // Edges adjacent to node 2
87 : {2, 3, 7}, // Edges adjacent to node 3
88 : {4, 99, 99}, // Edges adjacent to node 4
89 : {5, 99, 99}, // Edges adjacent to node 5
90 : {6, 99, 99}, // Edges adjacent to node 6
91 : {7, 99, 99} // Edges adjacent to node 7
92 : };
93 :
94 : // ------------------------------------------------------------
95 : // InfHex class member functions
96 0 : dof_id_type InfHex::key (const unsigned int s) const
97 : {
98 0 : libmesh_assert_less (s, this->n_sides());
99 :
100 : // The order of the node ids does not matter, they are sorted by the
101 : // compute_key() function.
102 0 : return this->compute_key(this->node_id(InfHex8::side_nodes_map[s][0]),
103 0 : this->node_id(InfHex8::side_nodes_map[s][1]),
104 0 : this->node_id(InfHex8::side_nodes_map[s][2]),
105 0 : this->node_id(InfHex8::side_nodes_map[s][3]));
106 : }
107 :
108 :
109 :
110 30845 : dof_id_type InfHex::low_order_key (const unsigned int s) const
111 : {
112 12250 : libmesh_assert_less (s, this->n_sides());
113 :
114 : // The order of the node ids does not matter, they are sorted by the
115 : // compute_key() function.
116 43095 : return this->compute_key(this->node_id(InfHex8::side_nodes_map[s][0]),
117 30845 : this->node_id(InfHex8::side_nodes_map[s][1]),
118 30845 : this->node_id(InfHex8::side_nodes_map[s][2]),
119 43095 : this->node_id(InfHex8::side_nodes_map[s][3]));
120 : }
121 :
122 :
123 :
124 0 : unsigned int InfHex::local_side_node(unsigned int side,
125 : unsigned int side_node) const
126 : {
127 0 : libmesh_assert_less (side, this->n_sides());
128 0 : libmesh_assert_less (side_node, InfHex8::nodes_per_side);
129 :
130 0 : return InfHex8::side_nodes_map[side][side_node];
131 : }
132 :
133 :
134 :
135 480 : unsigned int InfHex::local_edge_node(unsigned int edge,
136 : unsigned int edge_node) const
137 : {
138 160 : libmesh_assert_less (edge, this->n_edges());
139 160 : libmesh_assert_less (edge_node, InfHex8::nodes_per_edge);
140 :
141 480 : return InfHex8::edge_nodes_map[edge][edge_node];
142 : }
143 :
144 :
145 :
146 16296 : std::unique_ptr<Elem> InfHex::side_ptr (const unsigned int i)
147 : {
148 6484 : libmesh_assert_less (i, this->n_sides());
149 :
150 : // Return value
151 16296 : std::unique_ptr<Elem> face;
152 :
153 : // Think of a unit cube: (-1,1) x (-1,1) x (-1,1),
154 : // with (in general) the normals pointing outwards
155 16296 : switch (i)
156 : {
157 : // the face at z = -1
158 : // the base, where the infinite element couples to conventional
159 : // elements
160 5439 : case 0:
161 : {
162 : // Oops, here we are, claiming the normal of the face
163 : // elements point outwards -- and this is the exception:
164 : // For the side built from the base face,
165 : // the normal is pointing _into_ the element!
166 : // Why is that? - In agreement with build_side_ptr(),
167 : // which in turn _has_ to build the face in this
168 : // way as to enable the cool way \p InfFE re-uses \p FE.
169 5439 : face = std::make_unique<Quad4>();
170 5439 : break;
171 : }
172 :
173 : // These faces connect to other infinite elements.
174 10857 : case 1: // the face at y = -1
175 : case 2: // the face at x = 1
176 : case 3: // the face at y = 1
177 : case 4: // the face at x = -1
178 : {
179 10857 : face = std::make_unique<InfQuad4>();
180 10857 : break;
181 : }
182 :
183 0 : default:
184 0 : libmesh_error_msg("Invalid side i = " << i);
185 : }
186 :
187 : // Set the nodes
188 81480 : for (auto n : face->node_index_range())
189 91120 : face->set_node(n, this->node_ptr(InfHex8::side_nodes_map[i][n]));
190 :
191 16296 : return face;
192 0 : }
193 :
194 :
195 :
196 28689 : void InfHex::side_ptr (std::unique_ptr<Elem> & side,
197 : const unsigned int i)
198 : {
199 11399 : libmesh_assert_less (i, this->n_sides());
200 :
201 : // Think of a unit cube: (-1,1) x (-1,1) x (1,1)
202 28689 : switch (i)
203 : {
204 : // the base face
205 2203 : case 0:
206 : {
207 5545 : if (!side.get() || side->type() != QUAD4)
208 : {
209 6532 : side = this->side_ptr(i);
210 5424 : return;
211 : }
212 45 : break;
213 : }
214 :
215 : // connecting to another infinite element
216 9196 : case 1:
217 : case 2:
218 : case 3:
219 : case 4:
220 : {
221 23144 : if (!side.get() || side->type() != INFQUAD4)
222 : {
223 13002 : side = this->side_ptr(i);
224 10797 : return;
225 : }
226 4900 : break;
227 : }
228 :
229 0 : default:
230 0 : libmesh_error_msg("Invalid side i = " << i);
231 : }
232 :
233 17413 : side->subdomain_id() = this->subdomain_id();
234 :
235 : // Set the nodes
236 62340 : for (auto n : side->node_index_range())
237 69652 : side->set_node(n, this->node_ptr(InfHex8::side_nodes_map[i][n]));
238 : }
239 :
240 :
241 :
242 14268 : bool InfHex::is_child_on_side(const unsigned int c,
243 : const unsigned int s) const
244 : {
245 13176 : libmesh_assert_less (c, this->n_children());
246 13176 : libmesh_assert_less (s, this->n_sides());
247 :
248 14268 : return (s == 0 || c+1 == s || c == s%4);
249 : }
250 :
251 :
252 :
253 98172 : bool InfHex::is_edge_on_side (const unsigned int e,
254 : const unsigned int s) const
255 : {
256 98076 : libmesh_assert_less (e, this->n_edges());
257 98076 : libmesh_assert_less (s, this->n_sides());
258 :
259 98172 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
260 : }
261 :
262 :
263 :
264 288 : std::vector<unsigned int> InfHex::sides_on_edge(const unsigned int e) const
265 : {
266 96 : libmesh_assert_less(e, this->n_edges());
267 :
268 288 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
269 : }
270 :
271 :
272 : bool
273 24 : InfHex::is_flipped() const
274 : {
275 15 : return (triple_product(this->point(1)-this->point(0),
276 15 : this->point(3)-this->point(0),
277 42 : this->point(4)-this->point(0)) < 0);
278 : }
279 :
280 : std::vector<unsigned int>
281 630 : InfHex::edges_adjacent_to_node(const unsigned int n) const
282 : {
283 210 : libmesh_assert_less(n, this->n_nodes());
284 :
285 : // For vertices, we use the InfHex::adjacent_edges_map with
286 : // appropriate "trimming" based on whether the vertices are "at
287 : // infinity" or not. Otherwise each of the mid-edge nodes is
288 : // adjacent only to the edge it is on, and the remaining nodes are
289 : // not adjacent to any edge.
290 630 : if (this->is_vertex(n))
291 : {
292 360 : auto trim = (n < 4) ? 0 : 2;
293 360 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n]) - trim};
294 : }
295 270 : else if (this->is_edge(n))
296 120 : return {n - this->n_vertices()};
297 :
298 50 : libmesh_assert(this->is_face(n) || this->is_internal(n));
299 100 : return {};
300 : }
301 :
302 63 : Real InfHex::quality (const ElemQuality q) const
303 : {
304 63 : switch (q)
305 : {
306 :
307 : /**
308 : * Compute the min/max diagonal ratio.
309 : * Source: CUBIT User's Manual.
310 : *
311 : * For infinite elements, we just only compute
312 : * the diagonal in the face...
313 : * Don't know whether this makes sense,
314 : * but should be a feasible way.
315 : */
316 0 : case DIAGONAL:
317 : {
318 : // Diagonal between node 0 and node 2
319 0 : const Real d02 = this->length(0,2);
320 :
321 : // Diagonal between node 1 and node 3
322 0 : const Real d13 = this->length(1,3);
323 :
324 : // Find the biggest and smallest diagonals
325 0 : const Real min = std::min(d02, d13);
326 0 : const Real max = std::max(d02, d13);
327 :
328 0 : libmesh_assert_not_equal_to (max, 0.0);
329 :
330 0 : return min / max;
331 :
332 : break;
333 : }
334 :
335 : /**
336 : * Minimum ratio of lengths derived from opposite edges.
337 : * Source: CUBIT User's Manual.
338 : *
339 : * For IFEMs, do this only for the base face...
340 : * Does this make sense?
341 : */
342 0 : case TAPER:
343 : {
344 :
345 : /**
346 : * Compute the side lengths.
347 : */
348 0 : const Real d01 = this->length(0,1);
349 0 : const Real d12 = this->length(1,2);
350 0 : const Real d23 = this->length(2,3);
351 0 : const Real d03 = this->length(0,3);
352 :
353 0 : std::vector<Real> edge_ratios(2);
354 :
355 : // Bottom
356 0 : edge_ratios[0] = std::min(d01, d23) / std::max(d01, d23);
357 0 : edge_ratios[1] = std::min(d03, d12) / std::max(d03, d12);
358 :
359 0 : return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
360 :
361 : break;
362 : }
363 :
364 :
365 : /**
366 : * Minimum edge length divided by max diagonal length.
367 : * Source: CUBIT User's Manual.
368 : *
369 : * And again, we mess around a bit, for the IFEMs...
370 : * Do this only for the base.
371 : */
372 0 : case STRETCH:
373 : {
374 : /**
375 : * Should this be a sqrt2, when we do this for the base only?
376 : */
377 0 : const Real sqrt3 = 1.73205080756888;
378 :
379 : /**
380 : * Compute the maximum diagonal in the base.
381 : */
382 0 : const Real d02 = this->length(0,2);
383 0 : const Real d13 = this->length(1,3);
384 0 : const Real max_diag = std::max(d02, d13);
385 :
386 0 : libmesh_assert_not_equal_to ( max_diag, 0.0 );
387 :
388 : /**
389 : * Compute the minimum edge length in the base.
390 : */
391 0 : std::vector<Real> edges(4);
392 0 : edges[0] = this->length(0,1);
393 0 : edges[1] = this->length(1,2);
394 0 : edges[2] = this->length(2,3);
395 0 : edges[3] = this->length(0,3);
396 :
397 0 : const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
398 0 : return sqrt3 * min_edge / max_diag ;
399 : }
400 :
401 :
402 : /**
403 : * I don't know what to do for this metric.
404 : * Maybe the base class knows...
405 : */
406 63 : default:
407 63 : return Elem::quality(q);
408 : }
409 : }
410 :
411 :
412 :
413 0 : std::pair<Real, Real> InfHex::qual_bounds (const ElemQuality) const
414 : {
415 0 : libmesh_not_implemented();
416 :
417 : // We'll never get here
418 : return {0., 0.};
419 : }
420 :
421 :
422 :
423 :
424 :
425 : const unsigned short int InfHex::_second_order_adjacent_vertices[8][2] =
426 : {
427 : { 0, 1}, // vertices adjacent to node 8
428 : { 1, 2}, // vertices adjacent to node 9
429 : { 2, 3}, // vertices adjacent to node 10
430 : { 0, 3}, // vertices adjacent to node 11
431 :
432 : { 4, 5}, // vertices adjacent to node 12
433 : { 5, 6}, // vertices adjacent to node 13
434 : { 6, 7}, // vertices adjacent to node 14
435 : { 4, 7} // vertices adjacent to node 15
436 : };
437 :
438 :
439 :
440 : const unsigned short int InfHex::_second_order_vertex_child_number[18] =
441 : {
442 : 99,99,99,99,99,99,99,99, // Vertices
443 : 0,1,2,0, // Edges
444 : 0,1,2,0,0, // Faces
445 : 0 // Interior
446 : };
447 :
448 :
449 :
450 : const unsigned short int InfHex::_second_order_vertex_child_index[18] =
451 : {
452 : 99,99,99,99,99,99,99,99, // Vertices
453 : 1,2,3,3, // Edges
454 : 5,6,7,7,2, // Faces
455 : 6 // Interior
456 : };
457 :
458 :
459 3875 : bool InfHex::contains_point (const Point & p, Real tol) const
460 : {
461 : // For infinite elements with linear base interpolation:
462 : // make use of the fact that infinite elements do not
463 : // live inside the envelope. Use a fast scheme to
464 : // check whether point \p p is inside or outside
465 : // our relevant part of the envelope. Note that
466 : // this is not exclusive: only when the distance is less,
467 : // we are safe. Otherwise, we cannot say anything. The
468 : // envelope may be non-spherical, the physical point may lie
469 : // inside the envelope, outside the envelope, or even inside
470 : // this infinite element. Therefore if this fails,
471 : // fall back to the inverse_map()
472 3875 : const Point my_origin (this->origin());
473 :
474 : // determine the minimal distance of the base from the origin
475 : // Use norm_sq() instead of norm(), it is faster
476 12 : Point pt0_o(this->point(0) - my_origin);
477 6 : Point pt1_o(this->point(1) - my_origin);
478 6 : Point pt2_o(this->point(2) - my_origin);
479 6 : Point pt3_o(this->point(3) - my_origin);
480 3875 : const Real tmp_min_distance_sq = std::min(pt0_o.norm_sq(),
481 3875 : std::min(pt1_o.norm_sq(),
482 7732 : std::min(pt2_o.norm_sq(),
483 7756 : pt3_o.norm_sq())));
484 :
485 : // For a coarse grid, it is important to account for the fact
486 : // that the sides are not spherical, thus the closest point
487 : // can be closer than all edges.
488 : // This is an estimator using Pythagoras:
489 : const Real min_distance_sq = tmp_min_distance_sq
490 7738 : - .5*std::max((point(0)-point(2)).norm_sq(),
491 6202 : (point(1)-point(3)).norm_sq());
492 :
493 : // work with 1% allowable deviation. We can still fall
494 : // back to the InfFE::inverse_map()
495 3875 : const Real conservative_p_dist_sq = 1.01 * (Point(p - my_origin).norm_sq());
496 :
497 :
498 :
499 3875 : if (conservative_p_dist_sq < min_distance_sq)
500 : {
501 : // the physical point is definitely not contained in the element
502 6 : return false;
503 : }
504 :
505 : // this captures the case that the point is not (almost) in the direction of the element.:
506 : // first, project the problem onto the unit sphere:
507 0 : Point p_o(p - my_origin);
508 3810 : pt0_o /= pt0_o.norm();
509 3810 : pt1_o /= pt1_o.norm();
510 3810 : pt2_o /= pt2_o.norm();
511 3810 : pt3_o /= pt3_o.norm();
512 3810 : p_o /= p_o.norm();
513 :
514 :
515 : // now, check if it is in the projected face; using that the diagonal contains
516 : // the largest distance between points in it
517 7620 : Real max_h = std::max((pt0_o - pt2_o).norm_sq(),
518 5606 : (pt1_o - pt3_o).norm_sq())*1.01;
519 :
520 2644 : if ((p_o - pt0_o).norm_sq() > max_h ||
521 2024 : (p_o - pt1_o).norm_sq() > max_h ||
522 5171 : (p_o - pt2_o).norm_sq() > max_h ||
523 0 : (p_o - pt3_o).norm_sq() > max_h )
524 : {
525 : // the physical point is definitely not contained in the element
526 2820 : return false;
527 : }
528 :
529 : // Declare a basic FEType. Will use default in the base,
530 : // and something else (not important) in radial direction.
531 990 : FEType fe_type(default_order());
532 :
533 990 : const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
534 0 : tol, false);
535 :
536 990 : return this->on_reference_element(mapped_point, tol);
537 : }
538 :
539 :
540 :
541 1009 : bool InfHex::on_reference_element(const Point & p,
542 : const Real eps) const
543 : {
544 6 : const Real & xi = p(0);
545 6 : const Real & eta = p(1);
546 6 : const Real & zeta = p(2);
547 :
548 : // The reference infhex is a [-1,1]^3.
549 1969 : return ((xi >= -1.-eps) &&
550 960 : (xi <= 1.+eps) &&
551 929 : (eta >= -1.-eps) &&
552 888 : (eta <= 1.+eps) &&
553 1903 : (zeta >= -1.-eps) &&
554 1015 : (zeta <= 1.+eps));
555 : }
556 :
557 :
558 :
559 : } // namespace libMesh
560 :
561 :
562 :
563 :
564 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
