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 : #include "libmesh/libmesh_config.h"
21 :
22 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 :
24 : // Local includes
25 : #include "libmesh/face_inf_quad4.h"
26 : #include "libmesh/fe_interface.h"
27 : #include "libmesh/fe_type.h"
28 : #include "libmesh/edge_edge2.h"
29 : #include "libmesh/edge_inf_edge2.h"
30 : #include "libmesh/enum_io_package.h"
31 : #include "libmesh/enum_order.h"
32 : #include "libmesh/inf_fe_map.h"
33 :
34 : namespace libMesh
35 : {
36 :
37 :
38 :
39 : // ------------------------------------------------------------
40 : // InfQuad4 class static member initializations
41 : const int InfQuad4::num_nodes;
42 : const int InfQuad4::nodes_per_side;
43 :
44 : const unsigned int InfQuad4::side_nodes_map[InfQuad4::num_sides][InfQuad4::nodes_per_side] =
45 : {
46 : {0, 1}, // Side 0
47 : {1, 3}, // Side 1
48 : {0, 2} // Side 2
49 : };
50 :
51 : #ifdef LIBMESH_ENABLE_AMR
52 :
53 : const Real InfQuad4::_embedding_matrix[InfQuad4::num_children][InfQuad4::num_nodes][InfQuad4::num_nodes] =
54 : {
55 : // embedding matrix for child 0
56 : {
57 : // 0 1 2 3rd parent node
58 : {1.0, 0.0, 0.0, 0.0}, // 0th child node
59 : {0.5, 0.5, 0.0, 0.0}, // 1
60 : {0.0, 0.0, 1.0, 0.0}, // 2
61 : {0.0, 0.0, 0.5, 0.5} // 3
62 : },
63 :
64 : // embedding matrix for child 1
65 : {
66 : // 0 1 2 3
67 : {0.5, 0.5, 0.0, 0.0}, // 0
68 : {0.0, 1.0, 0.0, 0.0}, // 1
69 : {0.0, 0.0, 0.5, 0.5}, // 2
70 : {0.0, 0.0, 0.0, 1.0} // 3
71 : }
72 : };
73 :
74 : #endif
75 :
76 :
77 : // ------------------------------------------------------------
78 : // InfQuad4 class member functions
79 :
80 174 : bool InfQuad4::is_node_on_side(const unsigned int n,
81 : const unsigned int s) const
82 : {
83 66 : libmesh_assert_less (s, n_sides());
84 66 : return std::find(std::begin(side_nodes_map[s]),
85 174 : std::end(side_nodes_map[s]),
86 174 : n) != std::end(side_nodes_map[s]);
87 : }
88 :
89 : std::vector<unsigned>
90 75 : InfQuad4::nodes_on_side(const unsigned int s) const
91 : {
92 27 : libmesh_assert_less(s, n_sides());
93 102 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
94 : }
95 :
96 : std::vector<unsigned>
97 51 : InfQuad4::nodes_on_edge(const unsigned int e) const
98 : {
99 51 : return nodes_on_side(e);
100 : }
101 :
102 0 : bool InfQuad4::contains_point (const Point & p, Real tol) const
103 : {
104 : /*
105 : * make use of the fact that infinite elements do not
106 : * live inside the envelope. Use a fast scheme to
107 : * check whether point \p p is inside or outside
108 : * our relevant part of the envelope. Note that
109 : * this is not exclusive: the point may be outside
110 : * the envelope, but contained in another infinite element.
111 : * Therefore, if the distance is greater, do fall back
112 : * to the scheme of using inverse_map().
113 : */
114 0 : const Point my_origin (this->origin());
115 :
116 : /*
117 : * determine the minimal distance of the base from the origin
118 : * use norm_sq() instead of norm(), it is slightly faster
119 : */
120 0 : const Real min_distance_sq = std::min((Point(this->point(0)-my_origin)).norm_sq(),
121 0 : (Point(this->point(1)-my_origin)).norm_sq());
122 :
123 : /*
124 : * work with 1% allowable deviation. Can still fall
125 : * back to the InfFE::inverse_map()
126 : */
127 0 : const Real conservative_p_dist_sq = 1.01 * (Point(p-my_origin).norm_sq());
128 :
129 0 : if (conservative_p_dist_sq < min_distance_sq)
130 : {
131 : /*
132 : * the physical point is definitely not contained
133 : * in the element, return false.
134 : */
135 0 : return false;
136 : }
137 : else
138 : {
139 : /*
140 : * cannot say anything, fall back to the inverse_map()
141 : *
142 : * Declare a basic FEType. Will use default in the base,
143 : * and something else (not important) in radial direction.
144 : */
145 0 : FEType fe_type(default_order());
146 :
147 0 : const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
148 0 : tol, false);
149 :
150 0 : return this->on_reference_element(mapped_point, tol);
151 : }
152 : }
153 :
154 :
155 :
156 :
157 438 : Order InfQuad4::default_order() const
158 : {
159 438 : return FIRST;
160 : }
161 :
162 :
163 :
164 657 : std::unique_ptr<Elem> InfQuad4::build_side_ptr (const unsigned int i)
165 : {
166 : // libmesh_assert_less (i, this->n_sides());
167 :
168 657 : std::unique_ptr<Elem> edge;
169 :
170 657 : switch (i)
171 : {
172 227 : case 0:
173 : {
174 227 : edge = std::make_unique<Edge2>();
175 227 : break;
176 : }
177 :
178 : // adjacent to another infinite element
179 430 : case 1:
180 : case 2:
181 : {
182 430 : edge = std::make_unique<InfEdge2>();
183 430 : break;
184 : }
185 :
186 0 : default:
187 0 : libmesh_error_msg("Invalid side i = " << i);
188 : }
189 :
190 : // Set the nodes
191 1971 : for (auto n : edge->node_index_range())
192 1758 : edge->set_node(n, this->node_ptr(InfQuad4::side_nodes_map[i][n]));
193 :
194 657 : edge->set_interior_parent(this);
195 435 : edge->inherit_data_from(*this);
196 :
197 657 : return edge;
198 0 : }
199 :
200 :
201 :
202 42 : void InfQuad4::build_side_ptr (std::unique_ptr<Elem> & side,
203 : const unsigned int i)
204 : {
205 42 : this->side_ptr(side, i);
206 42 : side->set_interior_parent(this);
207 27 : side->inherit_data_from(*this);
208 42 : }
209 :
210 :
211 :
212 0 : void InfQuad4::connectivity(const unsigned int libmesh_dbg_var(sf),
213 : const IOPackage iop,
214 : std::vector<dof_id_type> & conn) const
215 : {
216 0 : libmesh_assert_less (sf, this->n_sub_elem());
217 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
218 :
219 0 : conn.resize(4);
220 :
221 0 : switch (iop)
222 : {
223 0 : case TECPLOT:
224 : {
225 0 : conn[0] = this->node_id(0)+1;
226 0 : conn[1] = this->node_id(1)+1;
227 0 : conn[2] = this->node_id(3)+1;
228 0 : conn[3] = this->node_id(2)+1;
229 0 : return;
230 : }
231 0 : case VTK:
232 : {
233 0 : conn[0] = this->node_id(0);
234 0 : conn[1] = this->node_id(1);
235 0 : conn[2] = this->node_id(3);
236 0 : conn[3] = this->node_id(2);
237 0 : return;
238 : }
239 0 : default:
240 0 : libmesh_error_msg("Unsupported IO package " << iop);
241 : }
242 : }
243 :
244 :
245 : ElemType
246 27 : InfQuad4::side_type (const unsigned int s) const
247 : {
248 9 : libmesh_assert_less (s, 3);
249 27 : if (s == 0)
250 9 : return EDGE2;
251 6 : return INFEDGE2;
252 : }
253 :
254 :
255 5 : void InfQuad4::flip(BoundaryInfo * boundary_info)
256 : {
257 2 : libmesh_assert(boundary_info);
258 :
259 5 : swap2nodes(0,1);
260 5 : swap2nodes(2,3);
261 2 : swap2neighbors(1,2);
262 5 : swap2boundarysides(1,2,boundary_info);
263 5 : swap2boundaryedges(1,2,boundary_info);
264 5 : }
265 :
266 : } // namespace libMesh
267 :
268 :
269 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
