Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
20 : #include "libmesh/libmesh_config.h"
21 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 :
23 : // Local includes cont'd
24 : #include "libmesh/face_inf_quad.h"
25 : #include "libmesh/edge_edge2.h"
26 : #include "libmesh/edge_inf_edge2.h"
27 : #include "libmesh/face_inf_quad4.h"
28 : #include "libmesh/enum_elem_quality.h"
29 :
30 : namespace libMesh
31 : {
32 :
33 :
34 : // ------------------------------------------------------------
35 : // InfQuad class static member initializations
36 : const int InfQuad::num_sides;
37 : const int InfQuad::num_children;
38 :
39 : // Note: we can omit initialization of the third entry of each row because
40 : // static variables are automatically zero-initialized.
41 : const Real InfQuad::_master_points[6][3] =
42 : {
43 : {-1, 0},
44 : {1, 0},
45 : {-1, 1},
46 : {1, 1},
47 : {0, 0},
48 : {0, 1}
49 : };
50 :
51 : const unsigned int InfQuad::adjacent_sides_map[/*num_vertices*/4][/*max_adjacent_sides*/2] =
52 : {
53 : {0, 2}, // Sides adjacent to node 0
54 : {0, 1}, // Sides adjacent to node 1
55 : {2, 99}, // Sides adjacent to node 2
56 : {1, 99} // Sides adjacent to node 3
57 : };
58 :
59 : // ------------------------------------------------------------
60 : // InfQuad class member functions
61 0 : dof_id_type InfQuad::key (const unsigned int s) const
62 : {
63 0 : libmesh_assert_less (s, this->n_sides());
64 :
65 : // The order of the node ids does not matter, they are sorted by the
66 : // compute_key() function.
67 0 : return this->compute_key(this->node_id(InfQuad4::side_nodes_map[s][0]),
68 0 : this->node_id(InfQuad4::side_nodes_map[s][1]));
69 : }
70 :
71 :
72 :
73 810 : dof_id_type InfQuad::low_order_key (const unsigned int s) const
74 : {
75 324 : libmesh_assert_less (s, this->n_sides());
76 :
77 : // The order of the node ids does not matter, they are sorted by the
78 : // compute_key() function.
79 1134 : return this->compute_key(this->node_id(InfQuad4::side_nodes_map[s][0]),
80 1134 : this->node_id(InfQuad4::side_nodes_map[s][1]));
81 : }
82 :
83 :
84 :
85 114 : unsigned int InfQuad::local_side_node(unsigned int side,
86 : unsigned int side_node) const
87 : {
88 38 : libmesh_assert_less (side, this->n_sides());
89 38 : libmesh_assert_less (side_node, 2);
90 :
91 114 : return InfQuad4::side_nodes_map[side][side_node];
92 : }
93 :
94 :
95 :
96 228 : unsigned int InfQuad::local_edge_node(unsigned int edge,
97 : unsigned int edge_node) const
98 : {
99 228 : return local_side_node(edge, edge_node);
100 : }
101 :
102 :
103 :
104 112 : std::unique_ptr<Elem> InfQuad::side_ptr (const unsigned int i)
105 : {
106 44 : libmesh_assert_less (i, this->n_sides());
107 :
108 : // Return value
109 112 : std::unique_ptr<Elem> edge;
110 :
111 112 : switch (i)
112 : {
113 31 : case 0: // base face
114 : {
115 31 : edge = std::make_unique<Edge2>();
116 31 : break;
117 : }
118 :
119 81 : case 1: // adjacent to another infinite element
120 : case 2: // adjacent to another infinite element
121 : {
122 81 : edge = std::make_unique<InfEdge2>();
123 81 : break;
124 : }
125 :
126 0 : default:
127 0 : libmesh_error_msg("Invalid side i = " << i);
128 : }
129 :
130 : // Set the nodes
131 336 : for (auto n : edge->node_index_range())
132 312 : edge->set_node(n, this->node_ptr(InfQuad4::side_nodes_map[i][n]));
133 :
134 112 : return edge;
135 0 : }
136 :
137 :
138 :
139 312 : void InfQuad::side_ptr (std::unique_ptr<Elem> & side,
140 : const unsigned int i)
141 : {
142 123 : libmesh_assert_less (i, this->n_sides());
143 :
144 312 : switch (i)
145 : {
146 : // the base face
147 9 : case 0:
148 : {
149 24 : if (!side.get() || side->type() != EDGE2)
150 : {
151 26 : side = this->side_ptr(i);
152 21 : return;
153 : }
154 1 : break;
155 : }
156 :
157 : // connecting to another infinite element
158 114 : case 1:
159 : case 2:
160 : {
161 288 : if (!side.get() || side->type() != INFEDGE2)
162 : {
163 74 : side = this->side_ptr(i);
164 61 : return;
165 : }
166 90 : break;
167 : }
168 :
169 0 : default:
170 0 : libmesh_error_msg("Invalid side i = " << i);
171 : }
172 :
173 321 : side->subdomain_id() = this->subdomain_id();
174 :
175 : // Set the nodes
176 690 : for (auto n : side->node_index_range())
177 642 : side->set_node(n, this->node_ptr(InfQuad4::side_nodes_map[i][n]));
178 : }
179 :
180 :
181 336 : bool InfQuad::is_child_on_side(const unsigned int c,
182 : const unsigned int s) const
183 : {
184 188 : libmesh_assert_less (c, this->n_children());
185 188 : libmesh_assert_less (s, this->n_sides());
186 :
187 440 : return (s == 0 ||
188 468 : (c == 0 && s == 2) ||
189 390 : (c == 1 && s == 1));
190 : }
191 :
192 :
193 : bool
194 16 : InfQuad::is_flipped() const
195 : {
196 : return (
197 : #if LIBMESH_DIM > 2
198 : // Don't bother outside the XY plane
199 22 : this->point(0)(2) && this->point(1)(2) &&
200 38 : this->point(2)(2) && this->point(3)(2) &&
201 : #endif
202 22 : ((this->point(1)(0)-this->point(0)(0))*
203 16 : (this->point(2)(1)-this->point(0)(1)) >
204 22 : (this->point(2)(0)-this->point(0)(0))*
205 22 : (this->point(1)(1)-this->point(0)(1))));
206 : }
207 :
208 :
209 : std::vector<unsigned int>
210 150 : InfQuad::edges_adjacent_to_node(const unsigned int n) const
211 : {
212 50 : libmesh_assert_less(n, this->n_nodes());
213 :
214 : // For vertices, we use the adjacent_sides_map, otherwise node
215 : // 4 is on side 0 and node 5 is not any any side.
216 150 : if (this->is_vertex(n))
217 : {
218 80 : auto trim = (n < 2) ? 0 : 1;
219 120 : return {std::begin(adjacent_sides_map[n]), std::end(adjacent_sides_map[n]) - trim};
220 : }
221 30 : else if (n == 4)
222 15 : return {0};
223 : else
224 10 : return {};
225 : }
226 :
227 :
228 30 : Real InfQuad::quality (const ElemQuality q) const
229 : {
230 30 : return Elem::quality(q); // Not implemented
231 : }
232 :
233 :
234 :
235 :
236 0 : std::pair<Real, Real> InfQuad::qual_bounds (const ElemQuality q) const
237 : {
238 0 : std::pair<Real, Real> bounds;
239 :
240 0 : switch (q)
241 : {
242 :
243 0 : case ASPECT_RATIO:
244 0 : bounds.first = 1.;
245 0 : bounds.second = 4.;
246 0 : break;
247 :
248 0 : case SKEW:
249 0 : bounds.first = 0.;
250 0 : bounds.second = 0.5;
251 0 : break;
252 :
253 0 : case TAPER:
254 0 : bounds.first = 0.;
255 0 : bounds.second = 0.7;
256 0 : break;
257 :
258 0 : case WARP:
259 0 : bounds.first = 0.9;
260 0 : bounds.second = 1.;
261 0 : break;
262 :
263 0 : case STRETCH:
264 0 : bounds.first = 0.25;
265 0 : bounds.second = 1.;
266 0 : break;
267 :
268 0 : case MIN_ANGLE:
269 0 : bounds.first = 45.;
270 0 : bounds.second = 90.;
271 0 : break;
272 :
273 0 : case MAX_ANGLE:
274 0 : bounds.first = 90.;
275 0 : bounds.second = 135.;
276 0 : break;
277 :
278 0 : case CONDITION:
279 0 : bounds.first = 1.;
280 0 : bounds.second = 4.;
281 0 : break;
282 :
283 0 : case JACOBIAN:
284 : case SCALED_JACOBIAN:
285 0 : bounds.first = 0.5;
286 0 : bounds.second = 1.;
287 0 : break;
288 :
289 0 : case SHEAR:
290 : case SHAPE:
291 : case SIZE:
292 0 : bounds.first = 0.3;
293 0 : bounds.second = 1.;
294 0 : break;
295 :
296 0 : case DISTORTION:
297 0 : bounds.first = 0.6;
298 0 : bounds.second = 1.;
299 0 : break;
300 :
301 0 : default:
302 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
303 0 : bounds.first = -1;
304 0 : bounds.second = -1;
305 : }
306 :
307 0 : return bounds;
308 : }
309 :
310 :
311 :
312 0 : bool InfQuad::on_reference_element(const Point & p,
313 : const Real eps) const
314 : {
315 0 : const Real & xi = p(0);
316 0 : const Real & eta = p(1);
317 :
318 : // The reference infquad is [-1,1]^2.
319 0 : return ((xi >= -1.-eps) &&
320 0 : (xi <= 1.+eps) &&
321 0 : (eta >= -1.-eps) &&
322 0 : (eta <= 1.+eps));
323 : }
324 :
325 :
326 :
327 : } // namespace libMesh
328 :
329 :
330 :
331 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
