form_element.h
Go to the documentation of this file.
#ifndef _RHEO_FORM_ELEMENT_H
#define _RHEO_FORM_ELEMENT_H
//
// for checking the spectral convergence exp(-ak)
// on the reference element
//
// author: Pierre.Saramito@imag.fr
//
// date: 1 october 2017
//
#include "rheolef/quadrature.h"
#include "rheolef/basis.h"
namespace rheolef {
// -----------------------------------------------------------------------------
// mass
// -----------------------------------------------------------------------------
template<class Value>
typename scalar_traits<Value>::type
product(const Value& a, const Value& b) { return a*b; }
template<class T>
T
product(const point_basic<T>& a, const point_basic<T>& b) { return dot(a,b); }
template<class T, class Value, class Basis>
void
build_mass_generic (
const Basis& b,
reference_element hat_K,
Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
{
typedef reference_element::size_type size_type;
const T eps = 1e3*std::numeric_limits<T>::epsilon();
quadrature_option qopt;
qopt.set_order(2*b.degree());
quadrature<Float> quad (qopt);
size_type loc_ndof = b.ndof (hat_K);
Eigen::Matrix<Value,Eigen::Dynamic,1> phi (loc_ndof);
mass.resize (loc_ndof, loc_ndof);
for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_K),
last_q = quad.end(hat_K); iter_q != last_q; iter_q++) {
b.evaluate (hat_K, (*iter_q).x, phi);
for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
T coeff_w = (*iter_q).w*product(phi[loc_idof],phi[loc_jdof]);
if (fabs(coeff_w) < eps) continue;
mass.coeffRef(loc_idof,loc_jdof) += coeff_w;
}}
}
// TODO: some zeros could appear as sum upon quadrature (e.g. Pk_sherwin):
// => do a second pass to eliminate it
mass.makeCompressed();
}
template<class T, class Basis>
void
build_mass (
const Basis& b,
reference_element hat_K,
Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
{
switch (b.valued_tag()) {
case space_constant::scalar: build_mass_generic<T,T,Basis> (b, hat_K, mass); break;
case space_constant::vector: build_mass_generic<T,point_basic<T>,Basis> (b, hat_K, mass); break;
default: error_macro("tensor-like mass_element test: not yet");
}
}
// -----------------------------------------------------------------------------
// mass summed on element boundary, for "Pkd[sides]" tests
// -----------------------------------------------------------------------------
template<class T, class Basis>
void
build_mass_bdr (
const Basis& b,
reference_element tilde_K,
Eigen::SparseMatrix<T,Eigen::RowMajor>& mass_bdr)
{
typedef reference_element::size_type size_type;
const T eps = 1e3*std::numeric_limits<T>::epsilon();
quadrature_option qopt;
qopt.set_order(2*b.degree());
quadrature<Float> quad (qopt);
size_type loc_ndof = b.ndof (tilde_K);
Eigen::Matrix<T,Eigen::Dynamic,1> phi (loc_ndof);
mass_bdr.resize (loc_ndof, loc_ndof);
for (size_type loc_isid = 0, loc_nsid = tilde_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
side_information_type sid (tilde_K, loc_isid);
reference_element hat_S = sid.hat;
T tilde_Js = tilde_K.side_measure(loc_isid)/measure(hat_S);
for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_S),
last_q = quad.end(hat_S); iter_q != last_q; iter_q++) {
T tilde_wq = tilde_Js*(*iter_q).w;
if (b.option().is_restricted_to_sides()) {
b.evaluate_on_side (tilde_K, sid, (*iter_q).x, phi);
} else {
b.evaluate (tilde_K, (*iter_q).x, phi);
}
for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
T coeff_w = tilde_wq*product(phi[loc_idof],phi[loc_jdof]);
if (fabs(coeff_w) < eps) continue;
mass_bdr.coeffRef(loc_idof,loc_jdof) += coeff_w;
}}
}
}
}
// -----------------------------------------------------------------------------
// grad_grad
// -----------------------------------------------------------------------------
template<class T, class Basis>
void
build_grad_grad (
const Basis& b,
reference_element hat_K,
Eigen::SparseMatrix<T,Eigen::RowMajor>& grad_grad)
{
typedef reference_element::size_type size_type;
const T eps = 1e3*std::numeric_limits<T>::epsilon();
quadrature_option qopt;
qopt.set_order(2*b.degree());
quadrature<Float> quad (qopt);
size_type loc_ndof = b.ndof (hat_K);
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi (loc_ndof);
grad_grad.resize (loc_ndof, loc_ndof);
for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_K),
last_q = quad.end(hat_K); iter_q != last_q; iter_q++) {
b.grad_evaluate (hat_K, (*iter_q).x, grad_phi);
for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
T coeff_w = (*iter_q).w*dot(grad_phi[loc_idof],grad_phi[loc_jdof]);
if (fabs(coeff_w) < eps) continue;
grad_grad.coeffRef (loc_idof,loc_jdof) += coeff_w;
}}
}
}
} // namespace rheolef
#endif // _RHEO_FORM_ELEMENT_H
size_type
field::size_type size_type
Definition: branch.cc:425
phi
Float phi(const point &nu, Float a, Float b)
Definition: burgers_flux_godunov.icc:25
rheolef::point_basic
Definition: point.h:87
rheolef::quadrature
Definition: quadrature.h:186
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition: quadrature.h:195
rheolef::quadrature::end
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
rheolef::quadrature::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::side_measure
Float side_measure(size_type loc_isid) const
Definition: reference_element.cc:213
rheolef::reference_element::n_side
size_type n_side() const
Definition: reference_element.h:104
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::build_mass_bdr
void build_mass_bdr(const Basis &b, reference_element tilde_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass_bdr)
Definition: form_element.h:93
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
rheolef::measure
Float measure(reference_element hat_K)
Definition: reference_element.cc:65
rheolef::build_mass
void build_mass(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass)
Definition: form_element.h:77
rheolef::build_mass_generic
void build_mass_generic(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass)
Definition: form_element.h:48
rheolef::build_grad_grad
void build_grad_grad(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &grad_grad)
Definition: form_element.h:132
rheolef::product
scalar_traits< Value >::type product(const Value &a, const Value &b)
Definition: form_element.h:41
a
Definition: diffusion_isotropic.h:25
phi
Definition: phi.h:25
rheolef::scalar_traits::type
T type
Definition: point.h:324
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
epsilon
Float epsilon
Definition: transmission_error.cc:25