quadrature_rep.cc
Go to the documentation of this file.
#include "rheolef/quadrature.h"
#include "rheolef/reference_element_face_transformation.h"
#include "rheolef/rheostream.h"
#include <cctype> // for std::isdigit
namespace rheolef {
using namespace std;
// ----------------------------------------------------------------------------
// quadrature options
// ----------------------------------------------------------------------------
static const char* static_family_name[quadrature_option::max_family+1] = {
"gauss" ,
"gauss_lobatto" ,
"gauss_radau" ,
"middle_edge" ,
"superconvergent" ,
"equispaced" ,
"undefined"
};
static
quadrature_option::family_type
family_name2type (string name)
{
for (size_t i = 0; i < quadrature_option::max_family; ++i) {
if (static_family_name[i] == name) {
return (quadrature_option::family_type) (i);
}
}
error_macro ("unexpected quadrature family name `" << name << "'");
return quadrature_option::max_family;
}
void
quadrature_option::set_family (string name)
{
set_family (family_name2type (name));
}
string
quadrature_option::get_family_name() const
{
check_macro (_family >= 0 && _family <= max_family,
"unexpected quadrature family number = " << _family);
return static_family_name[_family];
}
string
quadrature_option::name() const
{
if (get_family() == default_family && get_order() == default_order) {
return "";
} else {
return get_family_name() + "(" + itos(get_order()) + ")";
}
#ifdef TO_CLEAN
return get_family_name() + "(" + itos(long(get_order())) + ")";
#endif // TO_CLEAN
}
void
quadrature_option::reset (const std::string& name)
{
if (name == "") {
// set options to default:
set_family (default_family);
set_order (default_order);
return;
}
size_type i = name.find ('(');
string family = name.substr (0, i);
size_type i1 = i+1, l1 = name.size()-i1-1;
check_macro (name.size() > i1+1,
"invalid quadrature name \""<<name<<"\"; expect e.g. \"gauss(3)\"");
string order_str = name.substr (i1, l1);
for (size_type i = 0; i < order_str.size(); ++i) {
// accepts also e.g. "gauss(-1)" as default initialization
check_macro (std::isdigit(order_str[i]) || order_str[i] == '-',
"invalid quadrature name \""<<name<<"\"; expect e.g. \"gauss(3)\"");
}
size_type order = atoi(order_str.c_str());
set_family (family);
set_order (order);
}
// ----------------------------------------------------------------------------
// quadrature on a specific reference element
// ----------------------------------------------------------------------------
template<class T>
void
quadrature_on_geo<T>::initialize (reference_element hat_K, quadrature_option opt)
{
base::resize (0);
switch (hat_K.variant()) {
case reference_element::p: init_point (opt); break;
case reference_element::e: init_edge (opt); break;
case reference_element::t: init_triangle (opt); break;
case reference_element::q: init_square (opt); break;
case reference_element::T: init_tetrahedron (opt); break;
case reference_element::P: init_prism (opt); break;
case reference_element::H: init_hexahedron (opt); break;
default:
fatal_macro ("quadrature formula on element type `"
<< hat_K.name() << "' are not yet implemented.");
}
}
template<class T>
void
quadrature_on_geo<T>::get_nodes (Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
{
size_type nq = base::size();
node.resize (nq);
for (size_type q = 0; q < nq; ++q) {
node[q] = base::operator[](q).x;
}
}
// ----------------------------------------------------------------------------
// quadrature on the full set of reference elements
// ----------------------------------------------------------------------------
template<class T>
quadrature_rep<T>::quadrature_rep (quadrature_option opt)
: _options(opt),
_quad(),
_initialized (reference_element::max_variant, false)
{
}
template<class T>
quadrature_rep<T>::quadrature_rep (const std::string& name)
: _options(name),
_quad(),
_initialized (reference_element::max_variant, false)
{
}
template<class T>
quadrature_rep<T>::quadrature_rep (const quadrature_rep<T>& q)
: _options (q._options),
_quad (q._quad),
_initialized (q._initialized)
{
}
template<class T>
const quadrature_rep<T>&
quadrature_rep<T>::operator= (const quadrature_rep<T>& q) {
_options = q._options;
_quad = q._quad;
_initialized = q._initialized;
return *this;
}
template<class T>
string
quadrature_rep<T>::name() const
{
return _options.name();
}
template<class T>
void
quadrature_rep<T>::_initialize (reference_element hat_K) const
{
reference_element::variant_type K_type = hat_K.variant();
_quad[K_type].initialize (hat_K, _options);
_initialized [K_type] = true;
}
template<class T>
void
quadrature_rep<T>::get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
{
if (!_initialized [hat_K.variant()]) _initialize (hat_K);
_quad [hat_K.variant()].get_nodes (node);
}
template<class T>
typename quadrature_rep<T>::const_iterator
quadrature_rep<T>::begin (reference_element hat_K) const
{
reference_element::variant_type K_type = hat_K.variant();
if (!_initialized [K_type]) _initialize (hat_K);
return _quad[K_type].begin();
}
template<class T>
typename quadrature_rep<T>::const_iterator
quadrature_rep<T>::end (reference_element hat_K) const
{
reference_element::variant_type K_type = hat_K.variant();
if (!_initialized [K_type]) _initialize (hat_K);
return _quad[K_type].end();
}
template<class T>
typename quadrature_rep<T>::size_type
quadrature_rep<T>::size (reference_element hat_K) const
{
reference_element::variant_type K_type = hat_K.variant();
if (!_initialized [K_type]) _initialize (hat_K);
return _quad[K_type].size();
}
template<class T>
const weighted_point<T>&
quadrature_rep<T>::operator() (reference_element hat_K, size_type q) const
{
reference_element::variant_type K_type = hat_K.variant();
if (!_initialized [K_type]) _initialize (hat_K);
return _quad[K_type][q];
}
template<class T>
ostream&
operator<< (ostream& out, const quadrature_on_geo<T>& x)
{
out << setprecision (numeric_limits<T>::digits10)
<< x.size() << endl;
for (size_t r = 0; r < x.size(); r++)
out << x[r].x << "\t" << x[r].w << endl;
return out;
}
template<class T>
ostream&
operator<< (ostream& out, const quadrature_rep<T>& x)
{
out << "quadrature" << endl
<< x._options.get_family_name() << " " << x._options.get_order() << endl;
for (size_t i = 0; i != size_t(reference_element::max_variant); i++) {
reference_element::variant_type K_type = reference_element::variant_type(i);
if (! x._initialized [K_type]) continue;
reference_element hat_K (K_type);
out << "reference_element " << hat_K.name() << endl
<< x._quad[K_type];
}
out << "end_quadrature" << endl;
return out;
}
// ----------------------------------------------------------------------------
// quadrature : pointer interface
// ----------------------------------------------------------------------------
template <class T>
quadrature_rep<T>::~quadrature_rep()
{
if (_options.name() != "") {
persistent_table<quadrature<T>>::unload (name());
}
}
template<class T>
void
quadrature<T>::reset (const std::string& name)
{
if (name == "") {
base::operator= (new_macro(quadrature_rep<T>()));
} else {
base::operator= (persistent_table<quadrature<T>>::load (name));
}
}
template<class T>
const quadrature_option&
quadrature<T>::get_options() const
{
return base::data().get_options();
}
template<class T>
quadrature<T>::quadrature (quadrature_option opt)
: base(),
persistent_table<quadrature<T>>()
{
reset (opt.name());
}
template<class T>
quadrature<T>::quadrature (const std::string& name)
: base(),
persistent_table<quadrature<T>>()
{
reset (name);
}
template<class T>
void
quadrature<T>::set_order (size_type order)
{
quadrature_option qopt = get_options();
qopt.set_order (order);
reset (qopt.name());
}
template<class T>
void
quadrature<T>::set_family (family_type ft)
{
quadrature_option qopt = get_options();
qopt.set_family (ft);
reset (qopt.name());
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class quadrature_on_geo<Float>; \
template class quadrature_rep<Float>; \
template class quadrature<T>; \
template ostream& operator<< (ostream&, const quadrature_on_geo<T>&); \
template ostream& operator<< (ostream&, const quadrature_rep<T>&);
_RHEOLEF_instanciation(Float)
} // namespace rheolef
rheolef::quadrature_rep::quadrature_rep
quadrature_rep(quadrature_option opt=quadrature_option())
Definition: quadrature_rep.cc:135
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::io::out
@ out
Definition: rheostream.h:167
rheolef::weighted_point
Definition: quadrature.h:71
rheolef::quadrature::reset
void reset(const std::string &name)
Definition: quadrature_rep.cc:254
rheolef::point_basic
Definition: point.h:87
rheolef::quadrature_on_geo
Definition: quadrature.h:82
rheolef::quadrature_rep
Definition: quadrature.h:128
rheolef::quadrature_rep::operator()
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition: quadrature_rep.cc:210
rheolef::details::reset
void reset(T &x)
Definition: keller_details.h:31
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::quadrature_rep::get_nodes
void get_nodes(reference_element hat_K, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition: quadrature_rep.cc:179
rheolef::quadrature::set_order
void set_order(size_type order)
Definition: quadrature_rep.cc:284
rheolef::quadrature_rep::operator=
const quadrature_rep & operator=(const quadrature_rep< T > &q)
Definition: quadrature_rep.cc:157
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::quadrature_on_geo::initialize
void initialize(reference_element hat_K, quadrature_option opt)
Definition: quadrature_rep.cc:105
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::quadrature_on_geo::get_nodes
void get_nodes(Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition: quadrature_rep.cc:123
rheolef::quadrature::quadrature
quadrature(const std::string &name="")
Definition: quadrature_rep.cc:276
rheolef::quadrature::family_type
rep::family_type family_type
Definition: quadrature.h:194
rheolef::quadrature::get_options
const quadrature_option & get_options() const
Definition: quadrature_rep.cc:264
rheolef::persistent_table
see the persistent_table page for the full documentation
Definition: persistent_table.h:84
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::smart_pointer
see the smart_pointer page for the full documentation
Definition: smart_pointer.h:351
rheolef::quadrature::name
std::string name() const
Definition: quadrature.h:213
rheolef::reference_element::variant_type
size_type variant_type
Definition: reference_element.h:72
rheolef::quadrature_rep::_quad
std::array< quadrature_on_geo< T >, reference_element::max_variant > _quad
Definition: quadrature.h:175
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
rheolef::quadrature::set_family
void set_family(family_type ft)
Definition: quadrature_rep.cc:292
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition: quadrature.h:88
rheolef::quadrature_rep::size
size_type size(reference_element hat_K) const
Definition: quadrature_rep.cc:202
rheolef::quadrature_rep::_initialized
std::vector< bool > _initialized
Definition: quadrature.h:176
Float
see the Float page for the full documentation
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::quadrature_rep::_initialize
void _initialize(reference_element hat_K) const
Definition: quadrature_rep.cc:171
rheolef::quadrature_rep::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature_rep.cc:186
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::static_family_name
static const char * static_family_name[quadrature_option::max_family+1]
Definition: quadrature_rep.cc:31
load
void load(idiststream &in, Float &p, field &uh)
Definition: p_laplacian_post.cc:64
rheolef::quadrature_rep::_options
quadrature_option _options
Definition: quadrature.h:173
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::quadrature_rep::end
const_iterator end(reference_element hat_K) const
Definition: quadrature_rep.cc:194
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::quadrature
Definition: quadrature.h:185
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::quadrature_rep::name
std::string name() const
Definition: quadrature_rep.cc:165
rheolef::quadrature_rep::size_type
quadrature_on_geo< T >::size_type size_type
Definition: quadrature.h:133
rheolef::quadrature_rep::const_iterator
std::vector< weighted_point< T > >::const_iterator const_iterator
Definition: quadrature.h:135
rheolef::quadrature_rep::~quadrature_rep
~quadrature_rep()
Definition: quadrature_rep.cc:246
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::quadrature::size_type
rep::size_type size_type
Definition: quadrature.h:193