31 void hermite_compute_coefs(
const Elem * elem, std::vector<std::vector<Real>> & dxdxi
34 , std::vector<Real> & dydxi, std::vector<Real> & dzdeta, std::vector<Real> & dxdzeta,
35 std::vector<Real> & dzdxi, std::vector<Real> & dxdeta, std::vector<Real> & dydzeta
42 const int n_mapping_shape_functions =
48 for (
int p = 0; p != 2; ++p)
61 for (
int i = 0; i != n_mapping_shape_functions; ++i)
64 (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
66 (mapping_elem_type, mapping_order, i, 1, dofpt[p]);
68 (mapping_elem_type, mapping_order, i, 2, dofpt[p]);
73 dxdxi[0][p] += point_i(0) * ddxi;
74 dxdxi[1][p] += point_i(1) * ddeta;
75 dxdxi[2][p] += point_i(2) * ddzeta;
77 dydxi[p] += point_i(1) * ddxi;
78 dzdeta[p] += point_i(2) * ddeta;
79 dxdzeta[p] += point_i(0) * ddzeta;
80 dzdxi[p] += point_i(2) * ddxi;
81 dxdeta[p] += point_i(0) * ddeta;
82 dydzeta[p] += point_i(1) * ddzeta;
87 libmesh_assert(dxdxi[0][p]);
88 libmesh_assert(dxdxi[1][p]);
89 libmesh_assert(dxdxi[2][p]);
104 Real hermite_bases_3D (std::vector<unsigned int> & bases1D,
105 const std::vector<std::vector<Real>> & dxdxi,
113 unsigned int e = o-2;
149 libmesh_error_msg(
"Invalid basis node = " << i/8);
152 unsigned int basisnum = i%8;
159 coef = dxdxi[0][bases1D[0]];
160 bases1D[0] += 2;
break;
162 coef = dxdxi[1][bases1D[1]];
163 bases1D[1] += 2;
break;
165 coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
166 bases1D[0] += 2; bases1D[1] += 2;
break;
168 coef = dxdxi[2][bases1D[2]];
169 bases1D[2] += 2;
break;
171 coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
172 bases1D[0] += 2; bases1D[2] += 2;
break;
174 coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
175 bases1D[1] += 2; bases1D[2] += 2;
break;
177 coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
178 bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2;
break;
180 libmesh_error_msg(
"Invalid basisnum = " << basisnum);
184 else if (i < 64 + 12*4*e)
186 unsigned int basisnum = (i - 64) % (4*e);
187 switch ((i - 64) / (4*e))
190 bases1D[0] = basisnum / 4 + 4;
191 bases1D[1] = basisnum % 4 / 2 * 2;
192 bases1D[2] = basisnum % 2 * 2;
193 if (basisnum % 4 / 2)
199 bases1D[0] = basisnum % 4 / 2 * 2 + 1;
200 bases1D[1] = basisnum / 4 + 4;
201 bases1D[2] = basisnum % 2 * 2;
202 if (basisnum % 4 / 2)
208 bases1D[0] = basisnum / 4 + 4;
209 bases1D[1] = basisnum % 4 / 2 * 2 + 1;
210 bases1D[2] = basisnum % 2 * 2;
211 if (basisnum % 4 / 2)
217 bases1D[0] = basisnum % 4 / 2 * 2;
218 bases1D[1] = basisnum / 4 + 4;
219 bases1D[2] = basisnum % 2 * 2;
220 if (basisnum % 4 / 2)
226 bases1D[0] = basisnum % 4 / 2 * 2;
227 bases1D[1] = basisnum % 2 * 2;
228 bases1D[2] = basisnum / 4 + 4;
229 if (basisnum % 4 / 2)
235 bases1D[0] = basisnum % 4 / 2 * 2 + 1;
236 bases1D[1] = basisnum % 2 * 2;
237 bases1D[2] = basisnum / 4 + 4;
238 if (basisnum % 4 / 2)
244 bases1D[0] = basisnum % 4 / 2 * 2 + 1;
245 bases1D[1] = basisnum % 2 * 2 + 1;
246 bases1D[2] = basisnum / 4 + 4;
247 if (basisnum % 4 / 2)
253 bases1D[0] = basisnum % 4 / 2 * 2;
254 bases1D[1] = basisnum % 2 * 2 + 1;
255 bases1D[2] = basisnum / 4 + 4;
256 if (basisnum % 4 / 2)
262 bases1D[0] = basisnum / 4 + 4;
263 bases1D[1] = basisnum % 4 / 2 * 2;
264 bases1D[2] = basisnum % 2 * 2 + 1;
265 if (basisnum % 4 / 2)
271 bases1D[0] = basisnum % 4 / 2 * 2 + 1;
272 bases1D[1] = basisnum / 4 + 4;
273 bases1D[2] = basisnum % 2 * 2;
274 if (basisnum % 4 / 2)
280 bases1D[0] = basisnum / 4 + 4;
281 bases1D[1] = basisnum % 4 / 2 * 2 + 1;
282 bases1D[2] = basisnum % 2 * 2 + 1;
283 if (basisnum % 4 / 2)
289 bases1D[0] = basisnum % 4 / 2 * 2;
290 bases1D[1] = basisnum / 4 + 4;
291 bases1D[2] = basisnum % 2 * 2 + 1;
292 if (basisnum % 4 / 2)
298 libmesh_error_msg(
"Invalid basis node = " << (i - 64) / (4*e));
302 else if (i < 64 + 12*4*e + 6*2*e*e)
304 unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e);
305 switch ((i - 64 - 12*4*e) / (2*e*e))
310 bases1D[2] = basisnum % 2 * 2;
316 bases1D[1] = basisnum % 2 * 2;
322 bases1D[0] = basisnum % 2 * 2 + 1;
330 bases1D[1] = basisnum % 2 * 2 + 1;
336 bases1D[0] = basisnum % 2 * 2;
345 bases1D[2] = basisnum % 2 * 2 + 1;
350 libmesh_error_msg(
"Invalid basis node = " << (i - 64 - 12*4*e) / (2*e*e));
356 unsigned int basisnum = i - 64 - 12*4*e;
363 libmesh_assert(coef);
381 libmesh_error_msg(
"Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
390 const unsigned int i,
393 libmesh_assert(elem);
395 std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
398 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
399 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
402 hermite_compute_coefs(elem, dxdxi
404 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
423 libmesh_assert_less (i, 64);
425 std::vector<unsigned int> bases1D;
427 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
435 libmesh_error_msg(
"ERROR: Unsupported element type " << type);
440 libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
453 libmesh_error_msg(
"Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
462 const unsigned int i,
463 const unsigned int j,
466 libmesh_assert(elem);
467 libmesh_assert (j == 0 || j == 1 || j == 2);
469 std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
472 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
473 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
476 hermite_compute_coefs(elem, dxdxi
478 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
497 libmesh_assert_less (i, 64);
499 std::vector<unsigned int> bases1D;
501 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
524 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
529 libmesh_error_msg(
"ERROR: Unsupported element type " << type);
534 libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
543 const unsigned int i,
544 const unsigned int j,
547 libmesh_assert(elem);
549 std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
552 std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
553 std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
556 hermite_compute_coefs(elem, dxdxi
558 , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
577 libmesh_assert_less (i, 64);
579 std::vector<unsigned int> bases1D;
581 Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
622 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
627 libmesh_error_msg(
"ERROR: Unsupported element type " << type);
632 libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
The base class for all geometric element types.
const unsigned char cube_number_column[]
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
static const Real TOLERANCE
static Real hermite_raw_shape_second_deriv(const unsigned int basis_num, const Real xi)
static Real hermite_raw_shape(const unsigned int basis_num, const Real xi)
const unsigned char cube_number_row[]
const unsigned char square_number_column[]
const unsigned char cube_number_page[]
virtual unsigned int n_shape_functions() const override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned char square_number_row[]
static Real hermite_raw_shape_deriv(const unsigned int basis_num, const Real xi)
virtual Order default_order() const =0
virtual ElemType type() const =0
A geometric point in (x,y,z) space.
const Point & point(const unsigned int i) const
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)