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 : // 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 :
27 : namespace libMesh
28 : {
29 :
30 : // ------------------------------------------------------------
31 : // Hierarchic-specific implementations
32 :
33 : // Anonymous namespace for local helper functions
34 : namespace {
35 :
36 133616 : void side_hierarchic_nodal_soln(const Elem * elem,
37 : const Order /* order */,
38 : const std::vector<Number> & /* elem_soln */,
39 : std::vector<Number> & nodal_soln,
40 : const bool /*add_p_level*/)
41 : {
42 133616 : const unsigned int n_nodes = elem->n_nodes();
43 :
44 11156 : std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
45 133616 : nodal_soln.resize(n_nodes, 0);
46 :
47 : // We request nodal solutions when plotting, for consistency with
48 : // other elements; plotting 0 on non-sides makes sense in that
49 : // context.
50 : // libmesh_warning("Nodal solution requested for a side element; this makes no sense.");
51 133616 : } // side_hierarchic_nodal_soln()
52 :
53 :
54 7008 : void side_hierarchic_side_nodal_soln
55 : (const Elem * elem, const Order o,
56 : const unsigned int side,
57 : const std::vector<Number> & elem_soln,
58 : std::vector<Number> & nodal_soln_on_side,
59 : const bool /*add_p_level*/)
60 : {
61 : // Cheat here for now: perturb vertices toward the side center so as
62 : // to make the values well-defined.
63 : const std::vector<unsigned int> side_nodes =
64 7592 : elem->nodes_on_side(side);
65 1168 : const std::size_t n_side_nodes = side_nodes.size();
66 :
67 584 : const FEType fe_type{o, SIDE_HIERARCHIC};
68 :
69 584 : Point side_center;
70 53952 : for (auto n : side_nodes)
71 46944 : side_center += elem->master_point(n);
72 7008 : side_center /= n_side_nodes;
73 :
74 7008 : nodal_soln_on_side.resize(n_side_nodes);
75 53952 : for (auto i : make_range(n_side_nodes))
76 : {
77 46944 : const auto n = side_nodes[i];
78 46944 : Point master_p = elem->master_point(n);
79 3912 : master_p += TOLERANCE*TOLERANCE*(side_center-master_p);
80 :
81 : const unsigned int n_sf =
82 46944 : FEInterface::n_shape_functions(fe_type, elem);
83 :
84 46944 : nodal_soln_on_side[i] = 0;
85 1749024 : for (auto j : make_range(n_sf))
86 1843920 : nodal_soln_on_side[i] += elem_soln[j] *
87 1702080 : FEInterface::shape(fe_type, elem, j, master_p);
88 : }
89 7008 : }
90 :
91 :
92 :
93 34460958 : unsigned int side_hierarchic_n_dofs_at_node(const ElemType t,
94 : const Order o,
95 : const unsigned int n)
96 : {
97 34460958 : switch (t)
98 : {
99 28338 : case EDGE2:
100 : case EDGE3:
101 : case EDGE4:
102 28338 : if (n < 2)
103 1733 : return 1; // One per side
104 : else
105 9160 : return 0;
106 855812 : case QUAD8:
107 : case QUADSHELL8:
108 : case QUAD9:
109 : case QUADSHELL9:
110 855812 : if (n > 3 && n < 8)
111 368426 : return o+1;
112 : else
113 38375 : return 0;
114 2974714 : case HEX27:
115 2974714 : if (n > 19 && n < 26)
116 632127 : return (o+1)*(o+1); // (o+1)^2 per side
117 : else
118 152378 : return 0;
119 3765475 : case TRI6:
120 : case TRI7:
121 3765475 : if (n > 2 && n < 6)
122 1636634 : return o+1;
123 : else
124 176846 : return 0;
125 22267294 : case TET14:
126 22267294 : if (n > 9)
127 5939124 : return (o+1)*(o+2)/2;
128 : else
129 1176532 : return 0;
130 4569325 : case PRISM20:
131 : case PRISM21:
132 4569325 : if (n > 19)
133 6850 : return 0;
134 4461375 : if (n > 17)
135 416250 : return (o+1)*(o+2)/2;
136 4045125 : if (n > 14)
137 636075 : return (o+1)*(o+1);
138 204700 : return 0;
139 0 : case INVALID_ELEM:
140 0 : return 0;
141 : // Without side nodes on all sides we can't support side elements
142 0 : default:
143 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SIDE_HIERARCHIC FE family!");
144 : }
145 : } // side_hierarchic_n_dofs()
146 :
147 :
148 :
149 13194897 : unsigned int side_hierarchic_n_dofs_at_node(const Elem & e,
150 : const Order o,
151 : const unsigned int n)
152 : {
153 14069546 : return side_hierarchic_n_dofs_at_node(e.type(), o, n);
154 : }
155 :
156 :
157 :
158 6404559 : unsigned int side_hierarchic_n_dofs(const ElemType t, const Order o)
159 : {
160 6404559 : switch (t)
161 : {
162 653 : case EDGE2:
163 : case EDGE3:
164 : case EDGE4:
165 653 : return 2; // One per side
166 12294 : case QUAD8:
167 : case QUADSHELL8:
168 : case QUAD9:
169 : case QUADSHELL9:
170 91390 : return ((o+1)*4); // o+1 per side
171 12078 : case HEX27:
172 92795 : return ((o+1)*(o+1)*6); // (o+1)^2 per side
173 170092 : case TRI6:
174 : case TRI7:
175 1664588 : return ((o+1)*3); // o+1 per side
176 438178 : case TET14:
177 4457658 : return (o+1)*(o+2)*2; // 4 sides, each (o+1)(o+2)/2
178 15300 : case PRISM20:
179 : case PRISM21:
180 103375 : return (o+1)*(o+1)*3+(o+1)*(o+2); // 2 tris, (o+1)(o+2)/2; 3 quads
181 0 : case INVALID_ELEM:
182 0 : return 0;
183 : // Without side nodes on all sides we can't support side elements
184 0 : default:
185 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for HIERARCHIC FE family!");
186 : }
187 : } // side_hierarchic_n_dofs()
188 :
189 :
190 :
191 5882457 : unsigned int side_hierarchic_n_dofs(const Elem * e, const Order o)
192 : {
193 6404559 : return side_hierarchic_n_dofs(e->type(), o);
194 : }
195 :
196 :
197 : } // anonymous namespace
198 :
199 :
200 : // Instantiate nodal_soln() function for every dimension
201 133616 : LIBMESH_FE_NODAL_SOLN(SIDE_HIERARCHIC, side_hierarchic_nodal_soln)
202 :
203 : // side_nodal_soln() has to be manually defined here, since nodal_soln
204 : // isn't well-defined so we can't fall back on it.
205 : template <>
206 0 : void FE<0, SIDE_HIERARCHIC>::side_nodal_soln
207 : (const Elem *, const Order,
208 : const unsigned int,
209 : const std::vector<Number> &,
210 : std::vector<Number> &,
211 : bool,
212 : const unsigned)
213 : {
214 0 : libmesh_error_msg("No side variables in 0D!");
215 : }
216 :
217 : template <>
218 288 : void FE<1, SIDE_HIERARCHIC>::side_nodal_soln
219 : (const Elem *, const Order,
220 : const unsigned int side,
221 : const std::vector<Number> & elem_soln,
222 : std::vector<Number> & nodal_soln_on_side,
223 : const bool /*add_p_level*/,
224 : const unsigned)
225 : {
226 24 : libmesh_assert_less(side, 2);
227 288 : nodal_soln_on_side.resize(1);
228 312 : nodal_soln_on_side[0] = elem_soln[side];
229 288 : }
230 :
231 :
232 : template <>
233 2688 : void FE<2, SIDE_HIERARCHIC>::side_nodal_soln
234 : (const Elem * elem, const Order o,
235 : const unsigned int side,
236 : const std::vector<Number> & elem_soln,
237 : std::vector<Number> & nodal_soln_on_side,
238 : const bool add_p_level,
239 : const unsigned)
240 : {
241 224 : libmesh_assert_equal_to(elem->dim(), 2);
242 2688 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
243 : nodal_soln_on_side,
244 : add_p_level);
245 2688 : }
246 :
247 :
248 : template <>
249 4320 : void FE<3, SIDE_HIERARCHIC>::side_nodal_soln
250 : (const Elem * elem, const Order o,
251 : const unsigned int side,
252 : const std::vector<Number> & elem_soln,
253 : std::vector<Number> & nodal_soln_on_side,
254 : const bool add_p_level,
255 : const unsigned)
256 : {
257 360 : libmesh_assert_equal_to(elem->dim(), 3);
258 4320 : side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
259 : nodal_soln_on_side,
260 : add_p_level);
261 4320 : }
262 :
263 :
264 :
265 : // Full specialization of n_dofs() function for every dimension
266 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
267 0 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
268 0 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
269 0 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
270 :
271 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
272 1753 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
273 1755978 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
274 4646828 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
275 :
276 : // Full specialization of n_dofs_at_node() function for every dimension.
277 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); }
278 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); }
279 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); }
280 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); }
281 :
282 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); }
283 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); }
284 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); }
285 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); }
286 :
287 : // Full specialization of n_dofs_per_elem() function for every dimension.
288 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
289 0 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
290 0 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
291 0 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
292 :
293 0 : template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
294 8368 : template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
295 612432 : template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
296 1632064 : template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
297 :
298 : // Side FEMs are discontinuous from side to side
299 0 : template <> FEContinuity FE<0,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
300 11913 : template <> FEContinuity FE<1,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
301 132372 : template <> FEContinuity FE<2,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
302 388506 : template <> FEContinuity FE<3,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
303 :
304 : // Side Hierarchic FEMs are hierarchic (duh!)
305 0 : template <> bool FE<0,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
306 0 : template <> bool FE<1,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
307 0 : template <> bool FE<2,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
308 0 : template <> bool FE<3,SIDE_HIERARCHIC>::is_hierarchic() const { return true; }
309 :
310 : #ifdef LIBMESH_ENABLE_AMR
311 : // compute_constraints() specializations are only needed for 2 and 3D
312 : template <>
313 65168 : void FE<2,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
314 : DofMap & dof_map,
315 : const unsigned int variable_number,
316 : const Elem * elem)
317 65168 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
318 :
319 : template <>
320 145056 : void FE<3,SIDE_HIERARCHIC>::compute_constraints (DofConstraints & constraints,
321 : DofMap & dof_map,
322 : const unsigned int variable_number,
323 : const Elem * elem)
324 145056 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
325 : #endif // #ifdef LIBMESH_ENABLE_AMR
326 :
327 : // Hierarchic FEM shapes need reinit
328 0 : template <> bool FE<0,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
329 650 : template <> bool FE<1,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
330 1695850 : template <> bool FE<2,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
331 4466313 : template <> bool FE<3,SIDE_HIERARCHIC>::shapes_need_reinit() const { return true; }
332 :
333 : } // namespace libMesh
|
