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 :
20 : // Local includes
21 : #include "libmesh/elem.h"
22 : #include "libmesh/enum_to_string.h"
23 : #include "libmesh/fe.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/fe_macro.h"
26 : #include "libmesh/int_range.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 : // ------------------------------------------------------------
32 : // Hierarchic-specific implementations
33 :
34 : // Anonymous namespace for local helper functions
35 : namespace {
36 :
37 44372 : void side_hierarchic_nodal_soln(const Elem * elem,
38 : const Order order,
39 : const std::vector<Number> & elem_soln,
40 : std::vector<Number> & nodal_soln,
41 : const bool add_p_level)
42 : {
43 44372 : const unsigned int n_nodes = elem->n_nodes();
44 :
45 44372 : nodal_soln.assign(n_nodes, 0);
46 :
47 : // We request nodal solutions when plotting, for consistency with
48 : // other elements. Side solutions do not have unique values at element
49 : // nodes shared by multiple sides, so average the side values there.
50 : // libmesh_warning("Nodal solution requested for a side element; this makes no sense.");
51 55529 : std::vector<unsigned int> nodal_soln_count(n_nodes, 0);
52 11157 : const FEType fe_type{order, SIDE_HIERARCHIC};
53 22314 : std::vector<Number> nodal_soln_on_side;
54 :
55 205660 : for (const auto side : elem->side_index_range())
56 : {
57 : const std::vector<unsigned int> side_nodes =
58 201802 : elem->nodes_on_side(side);
59 :
60 161288 : FEInterface::side_nodal_soln(fe_type,
61 : elem,
62 : side,
63 : elem_soln,
64 : nodal_soln_on_side,
65 : add_p_level);
66 40514 : libmesh_assert_equal_to(nodal_soln_on_side.size(), side_nodes.size());
67 :
68 1058848 : for (const auto i : index_range(side_nodes))
69 : {
70 897560 : const auto n = side_nodes[i];
71 1122526 : nodal_soln[n] += nodal_soln_on_side[i];
72 1122526 : ++nodal_soln_count[n];
73 : }
74 : }
75 :
76 529412 : for (const auto n : index_range(nodal_soln))
77 606684 : if (nodal_soln_count[n])
78 477968 : nodal_soln[n] /= nodal_soln_count[n];
79 44372 : } // side_hierarchic_nodal_soln()
80 :
81 :
82 163496 : void side_hierarchic_side_nodal_soln
83 : (const Elem * elem, const Order o,
84 : const unsigned int side,
85 : const std::vector<Number> & elem_soln,
86 : std::vector<Number> & nodal_soln_on_side,
87 : const bool /*add_p_level*/)
88 : {
89 : // Cheat here for now: perturb vertices toward the side center so as
90 : // to make the values well-defined.
91 : const std::vector<unsigned int> side_nodes =
92 204562 : elem->nodes_on_side(side);
93 82132 : const std::size_t n_side_nodes = side_nodes.size();
94 :
95 41066 : const FEType fe_type{o, SIDE_HIERARCHIC};
96 :
97 41066 : Point side_center;
98 1076576 : for (auto n : side_nodes)
99 913080 : side_center += elem->master_point(n);
100 163496 : side_center /= n_side_nodes;
101 :
102 163496 : nodal_soln_on_side.resize(n_side_nodes);
103 1076576 : for (auto i : make_range(n_side_nodes))
104 : {
105 913080 : const auto n = side_nodes[i];
106 913080 : Point master_p = elem->master_point(n);
107 228846 : master_p += TOLERANCE*TOLERANCE*(side_center-master_p);
108 :
109 : const unsigned int n_sf =
110 913080 : FEInterface::n_shape_functions(fe_type, elem);
111 :
112 915384 : nodal_soln_on_side[i] = 0;
113 10254528 : for (auto j : make_range(n_sf))
114 11680266 : nodal_soln_on_side[i] += elem_soln[j] *
115 9341448 : FEInterface::shape(fe_type, elem, j, master_p);
116 : }
117 163496 : }
118 :
119 :
120 :
121 9237187 : unsigned int side_hierarchic_n_dofs_at_node(const ElemType t,
122 : const Order o,
123 : const unsigned int n)
124 : {
125 9237187 : switch (t)
126 : {
127 9069 : case EDGE2:
128 : case EDGE3:
129 : case EDGE4:
130 9069 : if (n < 2)
131 1733 : return 1; // One per side
132 : else
133 2946 : return 0;
134 247791 : case QUAD8:
135 : case QUADSHELL8:
136 : case QUAD9:
137 : case QUADSHELL9:
138 247791 : if (n > 3 && n < 8)
139 107044 : return o+1;
140 : else
141 38375 : return 0;
142 727541 : case HEX27:
143 727541 : if (n > 19 && n < 26)
144 157846 : return (o+1)*(o+1); // (o+1)^2 per side
145 : else
146 152378 : return 0;
147 1145801 : case TRI6:
148 : case TRI7:
149 1145801 : if (n > 2 && n < 6)
150 498254 : return o+1;
151 : else
152 176846 : return 0;
153 6062210 : case TET14:
154 6062210 : if (n > 9)
155 1648108 : return (o+1)*(o+2)/2;
156 : else
157 1176532 : return 0;
158 1044775 : case PRISM20:
159 : case PRISM21:
160 1044775 : if (n > 19)
161 6850 : return 0;
162 1018800 : if (n > 17)
163 97500 : return (o+1)*(o+2)/2;
164 921300 : if (n > 14)
165 150050 : return (o+1)*(o+1);
166 204700 : return 0;
167 0 : case INVALID_ELEM:
168 0 : return 0;
169 : // Without side nodes on all sides we can't support side elements
170 0 : default:
171 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SIDE_HIERARCHIC FE family!");
172 : }
173 : } // side_hierarchic_n_dofs()
174 :
175 :
176 :
177 2324139 : unsigned int side_hierarchic_n_dofs_at_node(const Elem & e,
178 : const Order o,
179 : const unsigned int n)
180 : {
181 3198788 : return side_hierarchic_n_dofs_at_node(e.type(), o, n);
182 : }
183 :
184 :
185 :
186 3105549 : unsigned int side_hierarchic_n_dofs(const ElemType t, const Order o)
187 : {
188 3105549 : switch (t)
189 : {
190 653 : case EDGE2:
191 : case EDGE3:
192 : case EDGE4:
193 653 : return 2; // One per side
194 18086 : case QUAD8:
195 : case QUADSHELL8:
196 : case QUAD9:
197 : case QUADSHELL9:
198 55468 : return ((o+1)*4); // o+1 per side
199 19850 : case HEX27:
200 64244 : return ((o+1)*(o+1)*6); // (o+1)^2 per side
201 209422 : case TRI6:
202 : case TRI7:
203 729438 : return ((o+1)*3); // o+1 per side
204 610210 : case TET14:
205 2219346 : return (o+1)*(o+2)*2; // 4 sides, each (o+1)(o+2)/2
206 15300 : case PRISM20:
207 : case PRISM21:
208 43100 : return (o+1)*(o+1)*3+(o+1)*(o+2); // 2 tris, (o+1)(o+2)/2; 3 quads
209 0 : case INVALID_ELEM:
210 0 : return 0;
211 : // Without side nodes on all sides we can't support side elements
212 0 : default:
213 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for HIERARCHIC FE family!");
214 : }
215 : } // side_hierarchic_n_dofs()
216 :
217 :
218 :
219 2358521 : unsigned int side_hierarchic_n_dofs(const Elem * e, const Order o)
220 : {
221 3105549 : return side_hierarchic_n_dofs(e->type(), o);
222 : }
223 :
224 :
225 : } // anonymous namespace
226 :
227 :
228 : // Instantiate nodal_soln() function for every dimension
229 44372 : LIBMESH_FE_NODAL_SOLN(SIDE_HIERARCHIC, side_hierarchic_nodal_soln)
230 :
231 : // side_nodal_soln() has to be manually defined here, since nodal_soln
232 : // isn't well-defined so we can't fall back on it.
233 : template <>
234 0 : void FE<0, SIDE_HIERARCHIC>::side_nodal_soln
235 : (const Elem *, const Order,
236 : const unsigned int,
237 : const std::vector<Number> &,
238 : std::vector<Number> &,
239 : bool,
240 : const unsigned)
241 : {
242 0 : libmesh_error_msg("No side variables in 0D!");
243 : }
244 :
245 : template <>
246 224 : void FE<1, SIDE_HIERARCHIC>::side_nodal_soln
247 : (const Elem *, const Order,
248 : const unsigned int side,
249 : const std::vector<Number> & elem_soln,
250 : std::vector<Number> & nodal_soln_on_side,
251 : const bool /*add_p_level*/,
252 : const unsigned)
253 : {
254 56 : libmesh_assert_less(side, 2);
255 224 : nodal_soln_on_side.resize(1);
256 280 : nodal_soln_on_side[0] = elem_soln[side];
257 224 : }
258 :
259 :
260 : template <>
261 60296 : void FE<2, SIDE_HIERARCHIC>::side_nodal_soln
262 : (const Elem * elem, const Order o,
263 : const unsigned int side,
264 : const std::vector<Number> & elem_soln,
265 : std::vector<Number> & nodal_soln_on_side,
266 : const bool add_p_level,
267 : const unsigned)
268 : {
269 15266 : libmesh_assert_equal_to(elem->dim(), 2);
270 60296 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
271 : nodal_soln_on_side,
272 : add_p_level);
273 60296 : }
274 :
275 :
276 : template <>
277 103200 : void FE<3, SIDE_HIERARCHIC>::side_nodal_soln
278 : (const Elem * elem, const Order o,
279 : const unsigned int side,
280 : const std::vector<Number> & elem_soln,
281 : std::vector<Number> & nodal_soln_on_side,
282 : const bool add_p_level,
283 : const unsigned)
284 : {
285 25800 : libmesh_assert_equal_to(elem->dim(), 3);
286 103200 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
287 : nodal_soln_on_side,
288 : add_p_level);
289 103200 : }
290 :
291 :
292 :
293 : // Full specialization of n_dofs() function for every dimension
294 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
295 0 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
296 0 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
297 0 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
298 :
299 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
300 953 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
301 784906 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
302 2319690 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
303 :
304 : // Full specialization of n_dofs_at_node() function for every dimension.
305 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
306 5868 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
307 960561 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
308 5071970 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
309 :
310 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
311 3201 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
312 433031 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
313 2762556 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
314 :
315 : // Full specialization of n_dofs_per_elem() function for every dimension.
316 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
317 0 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
318 0 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
319 0 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
320 :
321 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
322 2682 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
323 185704 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
324 438980 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
325 :
326 : // Side FEMs are discontinuous from side to side
327 0 : template <> FEContinuity FE<0,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
328 1155 : template <> FEContinuity FE<1,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
329 32836 : template <> FEContinuity FE<2,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
330 113926 : template <> FEContinuity FE<3,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
331 :
332 : // Side Hierarchic FEMs are hierarchic (duh!)
333 0 : template <> bool FE<0,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
334 0 : template <> bool FE<1,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
335 0 : template <> bool FE<2,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
336 0 : template <> bool FE<3,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
337 :
338 : #ifdef LIBMESH_ENABLE_AMR
339 : // compute_constraints() specializations are only needed for 2 and 3D
340 : template <>
341 21552 : void FE<2,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
342 : DofMap & dof_map,
343 : const unsigned int variable_number,
344 : const Elem * elem)
345 21552 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
346 :
347 : template <>
348 48352 : void FE<3,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
349 : DofMap & dof_map,
350 : const unsigned int variable_number,
351 : const Elem * elem)
352 48352 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
353 : #endif // #ifdef LIBMESH_ENABLE_AMR
354 :
355 : // Hierarchic FEM shapes need reinit
356 0 : template <> bool FE<0,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
357 225 : template <> bool FE<1,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
358 563382 : template <> bool FE<2,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
359 1486631 : template <> bool FE<3,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
360 :
361 : } // namespace libMesh
|
