Loading...
Searching...
No Matches
basis_fem_Pk_lagrange.cc
Go to the documentation of this file.
#include "basis_fem_Pk_lagrange.h"
#include "rheolef/rheostream.h"
#include "piola_fem_lagrange.h"
#include "equispaced.icc"
#include "warburton.icc"
#include "fekete.icc"
#include "Pk_get_local_idof_on_side.icc"
#include "basis_on_pointset_evaluate.icc"
#include "eigen_util.h"
namespace rheolef {
using namespace std;
// =========================================================================
// utilities
// =========================================================================
template <class T>
void
basis_fem_Pk_lagrange<T>::initialize_local_first (
size_type k,
bool is_continuous,
std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo_internal,
std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo,
std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo_internal,
std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo,
std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension_internal,
std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension,
std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension_internal,
std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension)
{
// 1) ndof_on_subgeo
if (k == size_t(-1)) {
// P_{-1} == empty polynomial space
for (size_type map_d = 0; map_d < 4; ++map_d) {
ndof_on_subgeo_internal [map_d].fill (0);
}
} else if (k == 0) {
// P0 initialization
for (size_type map_d = 0; map_d < 4; ++map_d) {
ndof_on_subgeo_internal [map_d].fill (0);
for (size_type variant = reference_element::first_variant_by_dimension (map_d);
variant < reference_element:: last_variant_by_dimension (map_d);
++variant) {
ndof_on_subgeo_internal [map_d] [variant] = 1;
}
}
} else { // k > 0
for (size_type map_d = 0; map_d < 4; ++map_d) {
reference_element::init_local_nnode_by_variant (k, ndof_on_subgeo_internal [map_d]);
// clean upper-dimensional subgeos, since init_local_nnode_by_variant do not consider dimension
for (size_type variant = reference_element::first_variant_by_dimension (map_d+1);
variant < reference_element::max_variant; ++variant) {
ndof_on_subgeo_internal [map_d] [variant] = 0;
}
}
}
// when discontinuous, fix, but conserve the subgeo structure in _internal:
base::_helper_make_discontinuous_ndof_on_subgeo (is_continuous, ndof_on_subgeo_internal, ndof_on_subgeo);
// 2) deduce automatically first_idof_by_dimension:
base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo_internal, first_idof_by_dimension_internal);
base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo, first_idof_by_dimension);
// 3) Lagrange nodes follow dofs:
nnod_on_subgeo_internal = ndof_on_subgeo_internal;
nnod_on_subgeo = ndof_on_subgeo;
first_inod_by_dimension_internal = first_idof_by_dimension_internal;
first_inod_by_dimension = first_idof_by_dimension;
}
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_Pk_lagrange<T>::~basis_fem_Pk_lagrange()
{
}
template<class T>
basis_fem_Pk_lagrange<T>::basis_fem_Pk_lagrange (size_type degree, const basis_option& sopt)
: basis_rep<T> (sopt),
_raw_basis(),
_hat_node(),
_vdm(),
_inv_vdm()
{
base::_name = base::standard_naming (family_name(), degree, base::_sopt);
string R = "?";
switch (base::_sopt.get_raw_polynomial()) {
case basis_option::monomial: R = "M"; break;
case basis_option::bernstein: R = "B"; break;
case basis_option::dubiner: R = "D"; break;
default: error_macro ("unsupported polynomial: "<<sopt.get_raw_polynomial_name());
}
_raw_basis = basis_raw_basic<T> (R+std::to_string(degree));
_initialize_cstor_sizes();
// piola FEM transformation:
typedef piola_fem_lagrange<T> piola_fem_type;
base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
}
template<class T>
bool
basis_fem_Pk_lagrange<T>::is_nodal() const
{
return true;
}
template<class T>
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
basis_fem_Pk_lagrange<T>::hat_node (reference_element hat_K) const
{
base::_initialize_data_guard (hat_K);
return _hat_node [hat_K.variant()];
}
template<class T>
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
basis_fem_Pk_lagrange<T>::vdm (reference_element hat_K) const
{
base::_initialize_data_guard (hat_K);
return _vdm [hat_K.variant()];
}
template<class T>
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
basis_fem_Pk_lagrange<T>::inv_vdm (reference_element hat_K) const
{
base::_initialize_data_guard (hat_K);
return _inv_vdm [hat_K.variant()];
}
template<class T>
void
basis_fem_Pk_lagrange<T>::_initialize_cstor_sizes() const
{
basis_fem_Pk_lagrange<T>::initialize_local_first (
degree(),
base::is_continuous(),
base::_ndof_on_subgeo_internal,
base::_ndof_on_subgeo,
base::_nnod_on_subgeo_internal,
base::_nnod_on_subgeo,
base::_first_idof_by_dimension_internal,
base::_first_idof_by_dimension,
base::_first_inod_by_dimension_internal,
base::_first_inod_by_dimension);
}
template<class T>
void
basis_fem_Pk_lagrange<T>::_initialize_data (reference_element hat_K) const
{
size_type k = degree();
size_type variant = hat_K.variant();
// nodes:
switch (base::_sopt.get_node()) {
case basis_option::equispaced:
pointset_lagrange_equispaced (hat_K, k, _hat_node[variant]); break;
case basis_option::warburton:
pointset_lagrange_warburton (hat_K, k, _hat_node[variant]); break;
case basis_option::fekete:
pointset_lagrange_fekete (hat_K, k, _hat_node[variant]); break;
default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
}
// vdm:
details::basis_on_pointset_evaluate (_raw_basis, hat_K, _hat_node[variant], _vdm[variant]);
check_macro (invert(_vdm[variant], _inv_vdm[variant]),
"unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
}
// evaluation of all basis functions at hat_x:
template<class T>
void
basis_fem_Pk_lagrange<T>::evaluate (
reference_element hat_K,
const point_basic<T>& hat_x,
Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
{
base::_initialize_data_guard (hat_K);
Eigen::Matrix<T,Eigen::Dynamic,1> raw_value;
_raw_basis.evaluate (hat_K, hat_x, raw_value);
value = _inv_vdm[hat_K.variant()].transpose()*raw_value;
}
// evaluate the gradient:
template<class T>
void
basis_fem_Pk_lagrange<T>::grad_evaluate (
reference_element hat_K,
const point_basic<T>& hat_x,
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
{
base::_initialize_data_guard (hat_K);
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm [hat_K.variant()];
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> raw_v_value;
_raw_basis.grad_evaluate (hat_K, hat_x, raw_v_value);
size_type loc_ndof = raw_v_value.size();
value.resize (loc_ndof);
// TODO: point-valued blas2 : value := trans(inv_vdm)*raw_value
for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
value [loc_idof] = point_basic<T>(0,0,0);
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
value [loc_idof] += inv_vdm (loc_jdof,loc_idof)*raw_v_value[loc_jdof];
}
}
}
// extract local dof-indexes on a side
template<class T>
typename basis_fem_Pk_lagrange<T>::size_type
basis_fem_Pk_lagrange<T>::local_ndof_on_side (
reference_element hat_K,
const side_information_type& sid) const
{
return base::ndof (sid.hat);
}
template<class T>
void
basis_fem_Pk_lagrange<T>::local_idof_on_side (
reference_element hat_K,
const side_information_type& sid,
Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
{
details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);
}
// dofs for a scalar-valued function
template<class T>
void
basis_fem_Pk_lagrange<T>::_compute_dofs (
reference_element hat_K,
const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
{
dof = f_xnod; // TODO: could avoid a physical copy: check it in compute_dof(f) with is_nodal()
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class basis_fem_Pk_lagrange<T>;
_RHEOLEF_instanciation(Float)
}// namespace rheolef
Pk_get_local_idof_on_side.icc
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
basis_fem_Pk_lagrange.h
basis_on_pointset_evaluate.icc
Float
see the Float page for the full documentation
rheolef::basis_fem_Pk_lagrange::~basis_fem_Pk_lagrange
~basis_fem_Pk_lagrange()
Definition basis_fem_Pk_lagrange.cc:95
rheolef::basis_fem_Pk_lagrange::_raw_basis
basis_raw_basic< T > _raw_basis
Definition basis_fem_Pk_lagrange.h:136
rheolef::basis_fem_Pk_lagrange::vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & vdm(reference_element hat_K) const
Definition basis_fem_Pk_lagrange.cc:136
rheolef::basis_fem_Pk_lagrange::size_type
reference_element::size_type size_type
Definition basis_fem_Pk_lagrange.h:58
rheolef::basis_fem_Pk_lagrange::local_idof_on_side
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) const
Definition basis_fem_Pk_lagrange.cc:232
rheolef::basis_fem_Pk_lagrange::degree
size_type degree() const
Definition basis_fem_Pk_lagrange.h:69
rheolef::basis_fem_Pk_lagrange::_initialize_cstor_sizes
void _initialize_cstor_sizes() const
Definition basis_fem_Pk_lagrange.cc:150
rheolef::basis_fem_Pk_lagrange::evaluate
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
Definition basis_fem_Pk_lagrange.cc:189
rheolef::basis_fem_Pk_lagrange::initialize_local_first
static void initialize_local_first(size_type k, bool is_continuous, std::array< std::array< size_type, reference_element::max_variant >, 4 > &ndof_on_subgeo_internal, std::array< std::array< size_type, reference_element::max_variant >, 4 > &ndof_on_subgeo, std::array< std::array< size_type, reference_element::max_variant >, 4 > &nnod_on_subgeo_internal, std::array< std::array< size_type, reference_element::max_variant >, 4 > &nnod_on_subgeo, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_idof_by_dimension_internal, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_idof_by_dimension, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_inod_by_dimension_internal, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_inod_by_dimension)
Definition basis_fem_Pk_lagrange.cc:41
rheolef::basis_fem_Pk_lagrange::local_ndof_on_side
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) const
Definition basis_fem_Pk_lagrange.cc:224
rheolef::basis_fem_Pk_lagrange::grad_evaluate
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value) const
Definition basis_fem_Pk_lagrange.cc:202
rheolef::basis_fem_Pk_lagrange::_compute_dofs
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
Definition basis_fem_Pk_lagrange.cc:242
rheolef::basis_fem_Pk_lagrange::basis_fem_Pk_lagrange
basis_fem_Pk_lagrange(size_type degree, const basis_option &sopt)
Definition basis_fem_Pk_lagrange.cc:99
rheolef::basis_fem_Pk_lagrange::hat_node
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
Definition basis_fem_Pk_lagrange.cc:129
rheolef::basis_fem_Pk_lagrange::family_name
virtual std::string family_name() const
Definition basis_fem_Pk_lagrange.h:68
rheolef::basis_fem_Pk_lagrange::is_nodal
bool is_nodal() const
Definition basis_fem_Pk_lagrange.cc:123
rheolef::basis_fem_Pk_lagrange::inv_vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & inv_vdm(reference_element hat_K) const
Definition basis_fem_Pk_lagrange.cc:143
rheolef::basis_fem_Pk_lagrange::_initialize_data
void _initialize_data(reference_element hat_K) const
Definition basis_fem_Pk_lagrange.cc:166
rheolef::basis_option
see the basis_option page for the full documentation
Definition basis_option.h:93
rheolef::basis_option::equispaced
@ equispaced
Definition basis_option.h:102
rheolef::basis_option::fekete
@ fekete
Definition basis_option.h:104
rheolef::basis_option::warburton
@ warburton
Definition basis_option.h:103
rheolef::basis_option::bernstein
@ bernstein
Definition basis_option.h:110
rheolef::basis_option::monomial
@ monomial
Definition basis_option.h:109
rheolef::basis_option::dubiner
@ dubiner
Definition basis_option.h:111
rheolef::basis_rep
Definition basis.h:209
rheolef::basis_rep::vector
vector
Definition basis.h:323
rheolef::piola_fem_lagrange
Definition piola_fem_lagrange.h:65
rheolef::point_basic
Definition point.h:86
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition reference_element.h:82
rheolef::reference_element::init_local_nnode_by_variant
static void init_local_nnode_by_variant(size_type order, std::array< size_type, reference_element::max_variant > &loc_nnod_by_variant)
Definition reference_element.cc:206
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition reference_element.h:150
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition reference_element.h:148
rheolef::reference_element::variant
variant_type variant() const
Definition reference_element.h:99
rheolef::reference_element::name
char name() const
Definition reference_element.h:100
error_macro
#define error_macro(message)
Definition dis_macros.h:49
eigen_util.h
equispaced.icc
fekete.icc
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details::basis_on_pointset_evaluate
void basis_on_pointset_evaluate(const Basis &b, const reference_element &hat_K, const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &vdm)
Definition basis_on_pointset_evaluate.icc:31
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::pointset_lagrange_warburton
void pointset_lagrange_warburton(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
Definition warburton.icc:574
rheolef::pointset_lagrange_fekete
void pointset_lagrange_fekete(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
Definition fekete.icc:42
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
rheolef::pointset_lagrange_equispaced
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
Definition equispaced.icc:44
std
STL namespace.
piola_fem_lagrange.h
rheolef::side_information_type
Definition reference_element_face_transformation.h:37
rheolef::side_information_type::hat
reference_element hat
Definition reference_element_face_transformation.h:43
warburton.icc