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 133628 : 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 133628 : const unsigned int n_nodes = elem->n_nodes();
44 :
45 133628 : 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 144785 : 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 619028 : for (const auto side : elem->side_index_range())
56 : {
57 : const std::vector<unsigned int> side_nodes =
58 525914 : elem->nodes_on_side(side);
59 :
60 485400 : 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 3182688 : for (const auto i : index_range(side_nodes))
69 : {
70 2697288 : const auto n = side_nodes[i];
71 2922254 : nodal_soln[n] += nodal_soln_on_side[i];
72 2922254 : ++nodal_soln_count[n];
73 : }
74 : }
75 :
76 1591820 : for (const auto n : index_range(nodal_soln))
77 1579836 : if (nodal_soln_count[n])
78 1436976 : nodal_soln[n] /= nodal_soln_count[n];
79 133628 : } // side_hierarchic_nodal_soln()
80 :
81 :
82 492024 : 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 533090 : 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 3235872 : for (auto n : side_nodes)
99 2743848 : side_center += elem->master_point(n);
100 492024 : side_center /= n_side_nodes;
101 :
102 492024 : nodal_soln_on_side.resize(n_side_nodes);
103 3235872 : for (auto i : make_range(n_side_nodes))
104 : {
105 2743848 : const auto n = side_nodes[i];
106 2743848 : 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 2743848 : FEInterface::n_shape_functions(fe_type, elem);
111 :
112 2746152 : nodal_soln_on_side[i] = 0;
113 30795840 : for (auto j : make_range(n_sf))
114 30390810 : nodal_soln_on_side[i] += elem_soln[j] *
115 28051992 : FEInterface::shape(fe_type, elem, j, master_p);
116 : }
117 492024 : }
118 :
119 :
120 :
121 34460958 : unsigned int side_hierarchic_n_dofs_at_node(const ElemType t,
122 : const Order o,
123 : const unsigned int n)
124 : {
125 34460958 : switch (t)
126 : {
127 28338 : case EDGE2:
128 : case EDGE3:
129 : case EDGE4:
130 28338 : if (n < 2)
131 1733 : return 1; // One per side
132 : else
133 9160 : return 0;
134 855812 : case QUAD8:
135 : case QUADSHELL8:
136 : case QUAD9:
137 : case QUADSHELL9:
138 855812 : if (n > 3 && n < 8)
139 368426 : return o+1;
140 : else
141 38375 : return 0;
142 2974714 : case HEX27:
143 2974714 : if (n > 19 && n < 26)
144 632127 : return (o+1)*(o+1); // (o+1)^2 per side
145 : else
146 152378 : return 0;
147 3765475 : case TRI6:
148 : case TRI7:
149 3765475 : if (n > 2 && n < 6)
150 1636634 : return o+1;
151 : else
152 176846 : return 0;
153 22267294 : case TET14:
154 22267294 : if (n > 9)
155 5939124 : return (o+1)*(o+2)/2;
156 : else
157 1176532 : return 0;
158 4569325 : case PRISM20:
159 : case PRISM21:
160 4569325 : if (n > 19)
161 6850 : return 0;
162 4461375 : if (n > 17)
163 416250 : return (o+1)*(o+2)/2;
164 4045125 : if (n > 14)
165 636075 : 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 13194897 : unsigned int side_hierarchic_n_dofs_at_node(const Elem & e,
178 : const Order o,
179 : const unsigned int n)
180 : {
181 14069546 : return side_hierarchic_n_dofs_at_node(e.type(), o, n);
182 : }
183 :
184 :
185 :
186 9101439 : unsigned int side_hierarchic_n_dofs(const ElemType t, const Order o)
187 : {
188 9101439 : 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 160930 : return ((o+1)*4); // o+1 per side
199 19850 : case HEX27:
200 186095 : return ((o+1)*(o+1)*6); // (o+1)^2 per side
201 209422 : case TRI6:
202 : case TRI7:
203 2134244 : return ((o+1)*3); // o+1 per side
204 610210 : case TET14:
205 6522042 : return (o+1)*(o+2)*2; // 4 sides, each (o+1)(o+2)/2
206 15300 : case PRISM20:
207 : case PRISM21:
208 103375 : 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 8354411 : unsigned int side_hierarchic_n_dofs(const Elem * e, const Order o)
220 : {
221 9101439 : 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 133628 : 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 672 : 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 672 : nodal_soln_on_side.resize(1);
256 728 : nodal_soln_on_side[0] = elem_soln[side];
257 672 : }
258 :
259 :
260 : template <>
261 182424 : 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 182424 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
271 : nodal_soln_on_side,
272 : add_p_level);
273 182424 : }
274 :
275 :
276 : template <>
277 309600 : 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 309600 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
287 : nodal_soln_on_side,
288 : add_p_level);
289 309600 : }
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 1753 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
301 2295174 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
302 6804512 : 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 18372 : 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 3160453 : 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 17212587 : 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 9966 : 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 1460834 : 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 12598746 : 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 8368 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
323 612432 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
324 1632064 : 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 11715 : template <> FEContinuity FE<1,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
329 130920 : template <> FEContinuity FE<2,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
330 387138 : 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 65168 : void FE<2,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
342 : DofMap & dof_map,
343 : const unsigned int variable_number,
344 : const Elem * elem)
345 65168 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
346 :
347 : template <>
348 145056 : void FE<3,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
349 : DofMap & dof_map,
350 : const unsigned int variable_number,
351 : const Elem * elem)
352 145056 : { 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 650 : template <> bool FE<1,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
358 1696558 : template <> bool FE<2,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
359 4467011 : template <> bool FE<3,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
360 :
361 : } // namespace libMesh
|
