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 : // 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 352 : bool InfQuad::is_child_on_side(const unsigned int c,
182 : const unsigned int s) const
183 : {
184 204 : libmesh_assert_less (c, this->n_children());
185 204 : libmesh_assert_less (s, this->n_sides());
186 :
187 352 : return (s == 0 || s == c+1);
188 : }
189 :
190 :
191 : bool
192 16 : InfQuad::is_flipped() const
193 : {
194 : return (
195 : #if LIBMESH_DIM > 2
196 : // Don't bother outside the XY plane
197 22 : this->point(0)(2) && this->point(1)(2) &&
198 38 : this->point(2)(2) && this->point(3)(2) &&
199 : #endif
200 22 : ((this->point(1)(0)-this->point(0)(0))*
201 16 : (this->point(2)(1)-this->point(0)(1)) >
202 22 : (this->point(2)(0)-this->point(0)(0))*
203 22 : (this->point(1)(1)-this->point(0)(1))));
204 : }
205 :
206 :
207 : std::vector<unsigned int>
208 150 : InfQuad::edges_adjacent_to_node(const unsigned int n) const
209 : {
210 50 : libmesh_assert_less(n, this->n_nodes());
211 :
212 : // For vertices, we use the adjacent_sides_map, otherwise node
213 : // 4 is on side 0 and node 5 is not any any side.
214 150 : if (this->is_vertex(n))
215 : {
216 80 : auto trim = (n < 2) ? 0 : 1;
217 120 : return {std::begin(adjacent_sides_map[n]), std::end(adjacent_sides_map[n]) - trim};
218 : }
219 30 : else if (n == 4)
220 15 : return {0};
221 : else
222 10 : return {};
223 : }
224 :
225 :
226 30 : Real InfQuad::quality (const ElemQuality q) const
227 : {
228 30 : return Elem::quality(q); // Not implemented
229 : }
230 :
231 :
232 :
233 :
234 0 : std::pair<Real, Real> InfQuad::qual_bounds (const ElemQuality q) const
235 : {
236 0 : std::pair<Real, Real> bounds;
237 :
238 0 : switch (q)
239 : {
240 :
241 0 : case ASPECT_RATIO:
242 0 : bounds.first = 1.;
243 0 : bounds.second = 4.;
244 0 : break;
245 :
246 0 : case SKEW:
247 0 : bounds.first = 0.;
248 0 : bounds.second = 0.5;
249 0 : break;
250 :
251 0 : case TAPER:
252 0 : bounds.first = 0.;
253 0 : bounds.second = 0.7;
254 0 : break;
255 :
256 0 : case WARP:
257 0 : bounds.first = 0.9;
258 0 : bounds.second = 1.;
259 0 : break;
260 :
261 0 : case STRETCH:
262 0 : bounds.first = 0.25;
263 0 : bounds.second = 1.;
264 0 : break;
265 :
266 0 : case MIN_ANGLE:
267 0 : bounds.first = 45.;
268 0 : bounds.second = 90.;
269 0 : break;
270 :
271 0 : case MAX_ANGLE:
272 0 : bounds.first = 90.;
273 0 : bounds.second = 135.;
274 0 : break;
275 :
276 0 : case CONDITION:
277 0 : bounds.first = 1.;
278 0 : bounds.second = 4.;
279 0 : break;
280 :
281 0 : case JACOBIAN:
282 : case SCALED_JACOBIAN:
283 0 : bounds.first = 0.5;
284 0 : bounds.second = 1.;
285 0 : break;
286 :
287 0 : case SHEAR:
288 : case SHAPE:
289 : case SIZE:
290 0 : bounds.first = 0.3;
291 0 : bounds.second = 1.;
292 0 : break;
293 :
294 0 : case DISTORTION:
295 0 : bounds.first = 0.6;
296 0 : bounds.second = 1.;
297 0 : break;
298 :
299 0 : default:
300 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
301 0 : bounds.first = -1;
302 0 : bounds.second = -1;
303 : }
304 :
305 0 : return bounds;
306 : }
307 :
308 :
309 :
310 0 : bool InfQuad::on_reference_element(const Point & p,
311 : const Real eps) const
312 : {
313 0 : const Real & xi = p(0);
314 0 : const Real & eta = p(1);
315 :
316 : // The reference infquad is [-1,1]^2.
317 0 : return ((xi >= -1.-eps) &&
318 0 : (xi <= 1.+eps) &&
319 0 : (eta >= -1.-eps) &&
320 0 : (eta <= 1.+eps));
321 : }
322 :
323 :
324 :
325 : } // namespace libMesh
326 :
327 :
328 :
329 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
