20 #include "libmesh/fe.h"    21 #include "libmesh/elem.h"    22 #include "libmesh/fe_interface.h"    23 #include "libmesh/number_lookups.h"    24 #include "libmesh/enum_to_string.h"    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
    44   const int n_mapping_shape_functions =
    52   for (
int p = 0; p != 2; ++p)
    65       for (
int i = 0; i != n_mapping_shape_functions; ++i)
    67           const Real ddxi = shape_deriv_ptr
    68             (map_fe_type, elem, i, 0, dofpt[p], 
false);
    69           const Real ddeta = shape_deriv_ptr
    70             (map_fe_type, elem, i, 1, dofpt[p], 
false);
    71           const Real ddzeta = shape_deriv_ptr
    72             (map_fe_type, elem, i, 2, dofpt[p], 
false);
    77           dxdxi[0][p] += point_i(0) * ddxi;
    78           dxdxi[1][p] += point_i(1) * ddeta;
    79           dxdxi[2][p] += point_i(2) * ddzeta;
    81           dydxi[p] += point_i(1) * ddxi;
    82           dzdeta[p] += point_i(2) * ddeta;
    83           dxdzeta[p] += point_i(0) * ddzeta;
    84           dzdxi[p] += point_i(2) * ddxi;
    85           dxdeta[p] += point_i(0) * ddeta;
    86           dydzeta[p] += point_i(1) * ddzeta;
    96       libmesh_assert_less (std::abs(dydxi[p]), 
TOLERANCE);
    97       libmesh_assert_less (std::abs(dzdeta[p]), 
TOLERANCE);
    98       libmesh_assert_less (std::abs(dxdzeta[p]), 
TOLERANCE);
    99       libmesh_assert_less (std::abs(dzdxi[p]), 
TOLERANCE);
   100       libmesh_assert_less (std::abs(dxdeta[p]), 
TOLERANCE);
   101       libmesh_assert_less (std::abs(dydzeta[p]), 
TOLERANCE);
   108 Real hermite_bases_3D (std::vector<unsigned int> & bases1D,
   109                        const std::vector<std::vector<Real>> & dxdxi,
   117   unsigned int e = o-2;
   153           libmesh_error_msg(
"Invalid basis node = " << i/8);
   156       unsigned int basisnum = i%8;
   163           coef = dxdxi[0][bases1D[0]];
   164           bases1D[0] += 2; 
break;
   166           coef = dxdxi[1][bases1D[1]];
   167           bases1D[1] += 2; 
break;
   169           coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
   170           bases1D[0] += 2; bases1D[1] += 2; 
break;
   172           coef = dxdxi[2][bases1D[2]];
   173           bases1D[2] += 2; 
break;
   175           coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
   176           bases1D[0] += 2; bases1D[2] += 2; 
break;
   178           coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
   179           bases1D[1] += 2; bases1D[2] += 2; 
break;
   181           coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
   182           bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2; 
break;
   184           libmesh_error_msg(
"Invalid basisnum = " << basisnum);
   188   else if (i < 64 + 12*4*e)
   190       unsigned int basisnum = (i - 64) % (4*e);
   191       switch ((i - 64) / (4*e))
   194           bases1D[0] = basisnum / 4 + 4;
   195           bases1D[1] = basisnum % 4 / 2 * 2;
   196           bases1D[2] = basisnum % 2 * 2;
   197           if (basisnum % 4 / 2)
   203           bases1D[0] = basisnum % 4 / 2 * 2 + 1;
   204           bases1D[1] = basisnum / 4 + 4;
   205           bases1D[2] = basisnum % 2 * 2;
   206           if (basisnum % 4 / 2)
   212           bases1D[0] = basisnum / 4 + 4;
   213           bases1D[1] = basisnum % 4 / 2 * 2 + 1;
   214           bases1D[2] = basisnum % 2 * 2;
   215           if (basisnum % 4 / 2)
   221           bases1D[0] = basisnum % 4 / 2 * 2;
   222           bases1D[1] = basisnum / 4 + 4;
   223           bases1D[2] = basisnum % 2 * 2;
   224           if (basisnum % 4 / 2)
   230           bases1D[0] = basisnum % 4 / 2 * 2;
   231           bases1D[1] = basisnum % 2 * 2;
   232           bases1D[2] = basisnum / 4 + 4;
   233           if (basisnum % 4 / 2)
   239           bases1D[0] = basisnum % 4 / 2 * 2 + 1;
   240           bases1D[1] = basisnum % 2 * 2;
   241           bases1D[2] = basisnum / 4 + 4;
   242           if (basisnum % 4 / 2)
   248           bases1D[0] = basisnum % 4 / 2 * 2 + 1;
   249           bases1D[1] = basisnum % 2 * 2 + 1;
   250           bases1D[2] = basisnum / 4 + 4;
   251           if (basisnum % 4 / 2)
   257           bases1D[0] = basisnum % 4 / 2 * 2;
   258           bases1D[1] = basisnum % 2 * 2 + 1;
   259           bases1D[2] = basisnum / 4 + 4;
   260           if (basisnum % 4 / 2)
   266           bases1D[0] = basisnum / 4 + 4;
   267           bases1D[1] = basisnum % 4 / 2 * 2;
   268           bases1D[2] = basisnum % 2 * 2 + 1;
   269           if (basisnum % 4 / 2)
   275           bases1D[0] = basisnum % 4 / 2 * 2 + 1;
   276           bases1D[1] = basisnum / 4 + 4;
   277           bases1D[2] = basisnum % 2 * 2;
   278           if (basisnum % 4 / 2)
   284           bases1D[0] = basisnum / 4 + 4;
   285           bases1D[1] = basisnum % 4 / 2 * 2 + 1;
   286           bases1D[2] = basisnum % 2 * 2 + 1;
   287           if (basisnum % 4 / 2)
   293           bases1D[0] = basisnum % 4 / 2 * 2;
   294           bases1D[1] = basisnum / 4 + 4;
   295           bases1D[2] = basisnum % 2 * 2 + 1;
   296           if (basisnum % 4 / 2)
   302           libmesh_error_msg(
"Invalid basis node = " << (i - 64) / (4*e));
   306   else if (i < 64 + 12*4*e + 6*2*e*e)
   308       unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e);
   309       switch ((i - 64 - 12*4*e) / (2*e*e))
   314           bases1D[2] = basisnum % 2 * 2;
   320           bases1D[1] = basisnum % 2 * 2;
   326           bases1D[0] = basisnum % 2 * 2 + 1;
   334           bases1D[1] = basisnum % 2 * 2 + 1;
   340           bases1D[0] = basisnum % 2 * 2;
   349           bases1D[2] = basisnum % 2 * 2 + 1;
   354           libmesh_error_msg(
"Invalid basis node = " << (i - 64 - 12*4*e) / (2*e*e));
   360       unsigned int basisnum = i - 64 - 12*4*e;
   385                           const unsigned int i,
   387                           const bool add_p_level)
   391   std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
   394   std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
   395   std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
   398   hermite_compute_coefs(elem, dxdxi
   400                         , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
   406   const Order totalorder =
   407     order + add_p_level*elem->
p_level();
   420               libmesh_assert_less (i, 64);
   422               std::vector<unsigned int> bases1D;
   424               Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
   437       libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
   448   libmesh_error_msg(
"Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
   455                           const unsigned int i,
   457                           const bool add_p_level)
   467                                 const unsigned int i,
   468                                 const unsigned int j,
   470                                 const bool add_p_level)
   475   std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
   478   std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
   479   std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
   482   hermite_compute_coefs(elem, dxdxi
   484                         , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
   490   const Order totalorder =
   491     order + add_p_level*elem->
p_level();
   504               libmesh_assert_less (i, 64);
   506               std::vector<unsigned int> bases1D;
   508               Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
   531                   libmesh_error_msg(
"Invalid shape function derivative j = " << j);
   541       libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
   553   libmesh_error_msg(
"Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
   561                                 const unsigned int i,
   562                                 const unsigned int j,
   564                                 const bool add_p_level)
   571 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES   577                                        const unsigned int i,
   578                                        const unsigned int j,
   580                                        const bool add_p_level)
   584   std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
   587   std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
   588   std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
   591   hermite_compute_coefs(elem, dxdxi
   593                         , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
   599   const Order totalorder =
   600     order + add_p_level*elem->
p_level();
   613               libmesh_assert_less (i, 64);
   615               std::vector<unsigned int> bases1D;
   617               Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
   658                   libmesh_error_msg(
"Invalid shape function derivative j = " << j);
   668       libmesh_error_msg(
"ERROR: Unsupported polynomial order " << totalorder);
   680   libmesh_error_msg(
"Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
   687                                        const unsigned int i,
   688                                        const unsigned int j,
   690                                        const bool add_p_level)
   695 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES Real(* shape_deriv_ptr)(const FEType fet, const Elem *elem, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Typedef for pointer to a function that returns FE shape function derivative values. 
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
ElemType
Defines an enum for geometric element types. 
Order
defines an enum for polynomial orders. 
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
static constexpr Real TOLERANCE
This is the base class from which all geometric element types are derived. 
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
OrderWrapper order
The approximation order of the element. 
The libMesh namespace provides an interface to certain functionality in the library. 
static Real hermite_raw_shape_second_deriv(const unsigned int basis_num, const Real xi)
1D hermite functions on unit interval 
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[]
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t, const ElemType t)
std::string enum_to_string(const T e)
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)
FEFamily
defines an enum for finite element families. 
virtual Order default_order() const =0
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real 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)
static FEFamily map_fe_type(const Elem &elem)