equispaced.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_EQUISPACED_ICC
#define _RHEOLEF_EQUISPACED_ICC
//
// compute the equispaced point set on the reference element
//
// author: Pierre.Saramito@imag.fr
//
// date: 2 september 2017
//
#include "rheolef/reference_element.h"
#include "reference_element_aux.icc"
#include "rheolef/compiler_eigen.h"
namespace rheolef {
// The optional argument "interior" specifies
// how many points from the boundary to omit.
// For example, on an "edge" reference element with order=2 and interior=0,
// this function will return the vertices and midpoint,
// but with interior=1, it will only return the midpoint.
template<class T>
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)
{
trace_macro("hat_K="<<hat_K.name()<<": order_in="<<order_in<<", internal="<<internal);
typedef reference_element::size_type size_type;
typedef point_basic<size_type> ilat;
if (order_in == 0) {
hat_xnod.resize(1);
reference_element_barycenter (hat_K, hat_xnod[0]);
trace_macro("barycenter: hat_xnod [0]="<<hat_xnod[0]);
return;
}
size_t d = hat_K.dimension();
// when only internal nodes: compute the effective "order" and the lattice size
size_t order = 0, size = 0;
switch (hat_K.variant()) {
case reference_element::p:
order = order_in;
size = reference_element::n_node(hat_K.variant(), order);
break;
case reference_element::e:
case reference_element::t:
case reference_element::T:
order = order_in + (d+1)*internal;
size = (order >= (d+1)*internal) ? reference_element::n_node(hat_K.variant(), order-(d+1)*internal) : 0;
break;
case reference_element::q:
case reference_element::H:
order = order_in + 2*internal;
size = (order >= 2*internal) ? reference_element::n_node(hat_K.variant(), order-2*internal) : 0;
break;
case reference_element::P:
check_macro (internal == 0, "internal: not yet");
order = order_in;
size = reference_element::n_node(hat_K.variant(), order);
break;
default:
error_macro ("unexpected element type `"<<hat_K.name()<<"'");
}
size_t first = internal;
size_t last = (order >= internal) ? order-internal : 0;
trace_macro("hat_K="<<hat_K.name()<<": order="<<order<<", first="<<first<<", last="<<last<<", size="<<size);
hat_xnod.resize (size);
switch (hat_K.variant()) {
case reference_element::p: {
hat_xnod.resize (1);
hat_xnod [0] = point_basic<T> (T(0));
break;
}
case reference_element::e: {
for (size_type i = first; i <= last; i++) {
size_type loc_idof = reference_element_e::ilat2loc_inod (order-(d+1)*internal, ilat(i-internal));
hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)));
trace_macro("e::lattice("<<i<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
}
break;
}
case reference_element::t: {
for (size_type j = first; j <= last; j++) {
for (size_type i = first; i+j <= last; i++) {
size_type loc_idof = reference_element_t::ilat2loc_inod (order-(d+1)*internal, ilat(i-internal,j-internal));
hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)));
trace_macro("t::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
}}
break;
}
case reference_element::q: {
for (size_type j = first; j <= last; j++) {
for (size_type i = first; i <= last; i++) {
size_type loc_idof = reference_element_q::ilat2loc_inod (order-2*internal, ilat(i-internal,j-internal));
hat_xnod [loc_idof] = point_basic<T> (-1+2*(T(int(i)))/T(int(order)), -1+2*T(int(j))/T(int(order)));
trace_macro("q::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
}}
break;
}
case reference_element::T: {
check_macro (internal == 0, "internal: not yet");
for (size_type k = first; k <= last; k++) {
for (size_type j = first; j+k <= last; j++) {
for (size_type i = first; i+j+k <= last; i++) {
size_type loc_idof = reference_element_T::ilat2loc_inod (order, ilat(i,j,k));
hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)), T(int(k))/T(int(order)));
}}}
break;
}
case reference_element::P: {
check_macro (internal == 0, "internal: not yet");
for (size_type k = first; k <= last; k++) {
for (size_type j = first; j <= last; j++) {
for (size_type i = first; i+j <= last; i++) {
size_type loc_idof = reference_element_P::ilat2loc_inod (order, ilat(i,j,k));
hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)), -1+2*T(int(k))/T(int(order)));
}}}
break;
}
case reference_element::H: {
check_macro (internal == 0, "internal: not yet");
for (size_type k = first; k <= last; k++) {
for (size_type j = first; j <= last; j++) {
for (size_type i = first; i <= last; i++) {
size_type loc_idof = reference_element_H::ilat2loc_inod (order, ilat(i,j,k));
hat_xnod [loc_idof] = point_basic<T> (-1+2*(T(int(i)))/T(int(order)), -1+2*T(int(j))/T(int(order)), -1+2*T(int(k))/T(int(order)));
}}}
break;
}
default: error_macro ("unexpected element type `"<<hat_K.name()<<"'");
}
}
} // namespace rheolef
#endif // _RHEOLEF_EQUISPACED_ICC
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element_P::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:695
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::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
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
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
reference_element_aux.icc
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:481
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
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::reference_element_H::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:790
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
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::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::reference_element_barycenter
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
Definition: reference_element_aux.icc:38
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71
T
Expr1::float_type T
Definition: field_expr.h:218
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111