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 : // C++ includes
20 :
21 : // Local includes
22 : #include "libmesh/fe.h"
23 : #include "libmesh/elem.h"
24 :
25 :
26 : namespace libMesh
27 : {
28 :
29 :
30 216 : LIBMESH_DEFAULT_VECTORIZED_FE(1,XYZ)
31 :
32 :
33 : template <>
34 16824 : Real FE<1,XYZ>::shape(const Elem * elem,
35 : const Order libmesh_dbg_var(order),
36 : const unsigned int i,
37 : const Point & point_in,
38 : const bool libmesh_dbg_var(add_p_level))
39 : {
40 1402 : libmesh_assert(elem);
41 1402 : libmesh_assert_less_equal (i, order + add_p_level * elem->p_level());
42 :
43 16824 : Point avg = elem->vertex_average();
44 16824 : Real max_distance = 0.;
45 65472 : for (const Point & p : elem->node_ref_range())
46 : {
47 52702 : const Real distance = std::abs(avg(0) - p(0));
48 48648 : max_distance = std::max(distance, max_distance);
49 : }
50 :
51 16824 : const Real x = point_in(0);
52 1402 : const Real xc = avg(0);
53 16824 : const Real dx = (x - xc)/max_distance;
54 :
55 : // monomials. since they are hierarchic we only need one case block.
56 16824 : switch (i)
57 : {
58 382 : case 0:
59 382 : return 1.;
60 :
61 4584 : case 1:
62 4584 : return dx;
63 :
64 3672 : case 2:
65 3672 : return dx*dx;
66 :
67 2604 : case 3:
68 2604 : return dx*dx*dx;
69 :
70 1380 : case 4:
71 1380 : return dx*dx*dx*dx;
72 :
73 0 : default:
74 0 : Real val = 1.;
75 0 : for (unsigned int index = 0; index != i; ++index)
76 0 : val *= dx;
77 0 : return val;
78 : }
79 : }
80 :
81 :
82 :
83 : template <>
84 0 : Real FE<1,XYZ>::shape(const ElemType,
85 : const Order,
86 : const unsigned int,
87 : const Point &)
88 : {
89 0 : libmesh_error_msg("XYZ polynomials require the element.");
90 : return 0.;
91 : }
92 :
93 :
94 :
95 :
96 : template <>
97 0 : Real FE<1,XYZ>::shape(const FEType fet,
98 : const Elem * elem,
99 : const unsigned int i,
100 : const Point & p,
101 : const bool add_p_level)
102 : {
103 0 : return FE<1,XYZ>::shape(elem, fet.order, i, p, add_p_level);
104 : }
105 :
106 :
107 : template <>
108 10344 : Real FE<1,XYZ>::shape_deriv(const Elem * elem,
109 : const Order libmesh_dbg_var(order),
110 : const unsigned int i,
111 : const unsigned int libmesh_dbg_var(j),
112 : const Point & point_in,
113 : const bool libmesh_dbg_var(add_p_level))
114 : {
115 862 : libmesh_assert(elem);
116 862 : libmesh_assert_less_equal (i, order + add_p_level * elem->p_level());
117 :
118 : // only d()/dxi in 1D!
119 :
120 862 : libmesh_assert_equal_to (j, 0);
121 :
122 10344 : Point avg = elem->vertex_average();
123 10344 : Real max_distance = 0.;
124 40032 : for (const Point & p : elem->node_ref_range())
125 : {
126 32162 : const Real distance = std::abs(avg(0) - p(0));
127 29688 : max_distance = std::max(distance, max_distance);
128 : }
129 :
130 10344 : const Real x = point_in(0);
131 862 : const Real xc = avg(0);
132 10344 : const Real dx = (x - xc)/max_distance;
133 :
134 : // monomials. since they are hierarchic we only need one case block.
135 10344 : switch (i)
136 : {
137 242 : case 0:
138 242 : return 0.;
139 :
140 2904 : case 1:
141 2904 : return 1./max_distance;
142 :
143 2232 : case 2:
144 2232 : return 2.*dx/max_distance;
145 :
146 1524 : case 3:
147 1524 : return 3.*dx*dx/max_distance;
148 :
149 780 : case 4:
150 780 : return 4.*dx*dx*dx/max_distance;
151 :
152 0 : default:
153 0 : Real val = i;
154 0 : for (unsigned int index = 1; index != i; ++index)
155 0 : val *= dx;
156 0 : return val/max_distance;
157 : }
158 : }
159 :
160 :
161 :
162 : template <>
163 0 : Real FE<1,XYZ>::shape_deriv(const ElemType,
164 : const Order,
165 : const unsigned int,
166 : const unsigned int,
167 : const Point &)
168 : {
169 0 : libmesh_error_msg("XYZ polynomials require the element.");
170 : return 0.;
171 : }
172 :
173 :
174 : template <>
175 0 : Real FE<1,XYZ>::shape_deriv(const FEType fet,
176 : const Elem * elem,
177 : const unsigned int i,
178 : const unsigned int j,
179 : const Point & p,
180 : const bool add_p_level)
181 : {
182 0 : return FE<1,XYZ>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
183 : }
184 :
185 :
186 :
187 :
188 :
189 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
190 :
191 :
192 :
193 : template <>
194 10344 : Real FE<1,XYZ>::shape_second_deriv(const Elem * elem,
195 : const Order libmesh_dbg_var(order),
196 : const unsigned int i,
197 : const unsigned int libmesh_dbg_var(j),
198 : const Point & point_in,
199 : const bool libmesh_dbg_var(add_p_level))
200 : {
201 862 : libmesh_assert(elem);
202 862 : libmesh_assert_less_equal (i, order + add_p_level * elem->p_level());
203 :
204 : // only d2()/dxi2 in 1D!
205 :
206 862 : libmesh_assert_equal_to (j, 0);
207 :
208 10344 : Point avg = elem->vertex_average();
209 10344 : Real max_distance = 0.;
210 40032 : for (const Point & p : elem->node_ref_range())
211 : {
212 32162 : const Real distance = std::abs(avg(0) - p(0));
213 29688 : max_distance = std::max(distance, max_distance);
214 : }
215 :
216 10344 : const Real x = point_in(0);
217 862 : const Real xc = avg(0);
218 10344 : const Real dx = (x - xc)/max_distance;
219 10344 : const Real dist2 = pow(max_distance,2.);
220 :
221 : // monomials. since they are hierarchic we only need one case block.
222 10344 : switch (i)
223 : {
224 484 : case 0:
225 : case 1:
226 484 : return 0.;
227 :
228 2232 : case 2:
229 2232 : return 2./dist2;
230 :
231 1524 : case 3:
232 1524 : return 6.*dx/dist2;
233 :
234 780 : case 4:
235 780 : return 12.*dx*dx/dist2;
236 :
237 0 : default:
238 0 : Real val = 2.;
239 0 : for (unsigned int index = 2; index != i; ++index)
240 0 : val *= (index+1) * dx;
241 0 : return val/dist2;
242 : }
243 : }
244 :
245 :
246 : template <>
247 0 : Real FE<1,XYZ>::shape_second_deriv(const ElemType,
248 : const Order,
249 : const unsigned int,
250 : const unsigned int,
251 : const Point &)
252 : {
253 0 : libmesh_error_msg("XYZ polynomials require the element.");
254 : return 0.;
255 : }
256 :
257 :
258 :
259 : template <>
260 0 : Real FE<1,XYZ>::shape_second_deriv(const FEType fet,
261 : const Elem * elem,
262 : const unsigned int i,
263 : const unsigned int j,
264 : const Point & p,
265 : const bool add_p_level)
266 : {
267 0 : return FE<1,XYZ>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
268 : }
269 :
270 : #endif
271 :
272 : } // namespace libMesh
|
