basis_fem_sides.cc
Go to the documentation of this file.
// There are four possibilities :
// 1) space Mh (omega["sides"], "P1d");
// => constructed side by side
// with local discontinuous connections between sides of an element
// along vertices (2d) or edges (3d)
// and global continuity for uni-valued approx along side
// status : supported in Rheolef
// used by HDG for "lambda"
// 2) space Mh (omega["sides"], "P1");
// => constructed side by side
// with local continuous connections between sides of an element
// and global continuity for uni-valued approx along side
// status : *NOT* supported in Rheolef
// 3) space Mhs (omega, "P1d[sides]");
// => constructed element by element, defined on each element boundary
// with local discontinuous connections between sides of an element
// along vertices (2d) or edges (3d)
// and global continuity for uni-valued approx along side
// status : supported in Rheolef
// used by HDG-POST for "lambda_s"
// 4) space Mhs (omega, "P1[sides,global_discontinuous]");
// => constructed element by element, defined on each element boundary
// with local continuous connections between sides of an element
// along vertices (2d) or edges (3d)
// and global discontinuity for uni-valued approx along side
// status : *NOT* supported in Rheolef
//
// The case 3 is the aim of this file: "Pkd[sides]"
//
#include "basis_fem_sides.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "equispaced.icc"
#include "warburton.icc"
#include "fekete.icc"
#include "eigen_util.h"
namespace rheolef {
using namespace std;
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_sides<T>::~basis_fem_sides()
{
}
template<class T>
basis_fem_sides<T>::basis_fem_sides (const basis_basic<T>& side_basis)
: basis_rep<T>(side_basis.option()),
_side_basis(side_basis),
_hat_node(),
_first_sid_inod(),
_first_sid_idof(),
_sid_value()
{
base::_sopt.set_restricted_to_sides (true);
base::_name = base::standard_naming (family_name(), family_index(), base::_sopt);
_initialize_cstor_sizes();
base::_piola_fem = _side_basis.get_piola_fem();
if (base::option().is_continuous()) {
// "Pkd[sices]" : discontinuous and multi-valued along sides connections (edges in 3d, vertives in 2d)
// "Pk[sides]" : could be continuous, but not yet supported
warning_macro("continuous approximation restricted to sides not yet supported");
error_macro ("HINT: replace \""<<base::name()<<"\" by \"P"<<degree()<<"d[sides]\"");
}
}
template<class T>
void
basis_fem_sides<T>::_initialize_cstor_sizes() const
{
size_type k = degree();
for (size_type d = 0; d < 4; ++d) {
base::_ndof_on_subgeo [d].fill (0);
base::_nnod_on_subgeo [d].fill (0);
if (d == 0) continue;
for (size_type subgeo_variant = reference_element::first_variant_by_dimension(d-1);
subgeo_variant < reference_element:: last_variant_by_dimension(d-1);
subgeo_variant++) {
// customized for _side_basis="Pkd" on subgeo-variant : all dofs in "_side_basis" are merged on (d-1) = dim(subgeo_variant)
// => all dofs are associated to the "d" dimension in the present element
// => cannot handle yet the case _side_basis="RTkd" here : it requires "_side_basis.ndof_on_subgeo_continuous(d,subgeo_variant)
base::_ndof_on_subgeo [d][subgeo_variant] = _side_basis.ndof_on_subgeo (d-1, subgeo_variant);
base::_nnod_on_subgeo [d][subgeo_variant] = _side_basis.nnod_on_subgeo (d-1, subgeo_variant);
}
}
bool is_continuous = false;
base::_helper_discontinuous_ndof_on_subgeo_inplace_change (is_continuous, base::_ndof_on_subgeo);
base::_helper_discontinuous_ndof_on_subgeo_inplace_change (is_continuous, base::_nnod_on_subgeo);
base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo, base::_first_idof_by_dimension);
base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo, base::_first_inod_by_dimension);
}
template<class T>
void
basis_fem_sides<T>::_initialize_data (reference_element tilde_K) const
{
size_type k = degree();
size_type d = tilde_K.dimension();
size_type variant = tilde_K.variant();
// nodes:
_first_sid_idof [variant].fill (0);
size_type loc_nnod = base::_first_inod_by_dimension [variant][d+1];
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& tilde_node = _hat_node [variant];
tilde_node.resize (loc_nnod);
size_type loc_first_sid_inod = 0;
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;
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& side_hat_node = _side_basis.hat_node (hat_S);
size_type loc_sid_nnod = side_hat_node.size();
for (size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
// piola transform side nodes from hat_S to partial tilde_K
size_type loc_inod = loc_first_sid_inod + loc_sid_inod;
tilde_node[loc_inod] = reference_element_face_transformation (tilde_K, sid, side_hat_node[loc_sid_inod]);
}
loc_first_sid_inod += loc_sid_nnod;
_first_sid_idof [variant][loc_isid+1] = loc_first_sid_inod;
}
_first_sid_inod [variant] = _first_sid_idof [variant];
}
template<class T>
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
basis_fem_sides<T>::hat_node (reference_element hat_K) const
{
base::_initialize_data_guard (hat_K);
return _hat_node [hat_K.variant()];
}
// evaluation on a side of all related basis functions at hat_x in hat_S
// NOTE: by convention, the value array has a full size for all sides dofs
// and is filled by zero on all others sides
template<class T>
void
basis_fem_sides<T>::evaluate_on_side (
reference_element tilde_K,
const side_information_type& sid,
const point_basic<T>& hat_x,
Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
{
base::_initialize_data_guard (tilde_K);
size_type variant = tilde_K.variant();
size_type loc_ndof = base::ndof (tilde_K);
value.resize (loc_ndof);
value.fill(T(0));
size_type loc_first_sid_idof = _first_sid_inod [variant][sid.loc_isid];
reference_element hat_S = sid.hat;
_side_basis.evaluate (hat_S, hat_x, _sid_value);
for (size_type loc_sid_idof = 0, loc_sid_ndof = _sid_value.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
size_type loc_idof = loc_first_sid_idof + loc_sid_idof;
value [loc_idof] = _sid_value [loc_sid_idof];
}
}
// extract local dof-indexes on a side
template<class T>
typename basis_fem_sides<T>::size_type
basis_fem_sides<T>::local_ndof_on_side (
reference_element tilde_K,
const side_information_type& sid) const
{
base::_initialize_data_guard (tilde_K);
return _first_sid_inod [tilde_K.variant()][sid.loc_isid+1]
- _first_sid_inod [tilde_K.variant()][sid.loc_isid];
}
template<class T>
void
basis_fem_sides<T>::local_idof_on_side (
reference_element tilde_K,
const side_information_type& sid,
Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
{
base::_initialize_data_guard (tilde_K);
size_type variant = tilde_K.variant();
size_type first_loc_inod = _first_sid_inod [variant][sid.loc_isid];
size_type loc_sid_nnod = _first_sid_inod [variant][sid.loc_isid+1]
- _first_sid_inod [variant][sid.loc_isid];
loc_idof.resize (loc_sid_nnod);
for (size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
loc_idof [loc_sid_inod] = first_loc_inod + loc_sid_inod;
}
}
// dofs for a scalar-valued function
template<class T>
void
basis_fem_sides<T>::_compute_dofs (
reference_element tilde_K,
const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
{
check_macro (is_nodal(), "_compute_dofs: only supported for nodal basis");
dof = f_xnod;
}
// -----------------------------------------------------------------------------
// visu: gnuplot in elevation for 2d (t,q)
// -----------------------------------------------------------------------------
template<class T>
void
basis_fem_sides<T>::put_scalar_valued (ostream& os, reference_element tilde_K) const
{
bool verbose = iorheo::getverbose(os);
bool clean = iorheo::getclean(os);
bool execute = iorheo::getexecute(os);
size_type nsub = iorheo::getsubdivide(os);
if (nsub <= 1) nsub = 20; // default value
size_type loc_ndof = base::ndof (tilde_K);
string basename = "basis-" + base::name() + "-" + tilde_K.name();
string filelist;
// --------------------
// .plot
// --------------------
string plot_name = basename + ".plot";
string gdat_name = basename + ".gdat";
ofstream plot (plot_name.c_str());
cerr << "! file \"" << plot_name << "\" created" << endl;
filelist += " " + plot_name;
size_type d = tilde_K.dimension();
check_macro (d==2, "unsupported dimension " << d);
plot << "gdat = \"" << gdat_name << "\"" << endl
<< "set colors classic" << endl
<< "set size square" << endl
<< "pause_duration = 1.5" << endl
<< "set view 55, 145" << endl
;
if (tilde_K.variant() == reference_element::t) {
plot << "set xrange [-0.1:1.1]" << endl
<< "set yrange [-0.1:1.1]" << endl
<< "set zrange [-1.1:1.1]" << endl
<< "set arrow from 0,0,0 to 1,0,0 nohead lc 0" << endl
<< "set arrow from 1,0,0 to 0,1,0 nohead lc 0" << endl
<< "set arrow from 0,1,0 to 0,0,0 nohead lc 0" << endl
;
} else {
plot << "set xrange [-1.1:1.1]" << endl
<< "set yrange [-1.1:1.1]" << endl
<< "set zrange [-1.1:1.1]" << endl
<< "set arrow from -1,-1,0 to 1,-1,0 nohead lc 0" << endl
<< "set arrow from 1,-1,0 to 1, 1,0 nohead lc 0" << endl
<< "set arrow from 1, 1,0 to -1, 1,0 nohead lc 0" << endl
<< "set arrow from -1, 1,0 to -1,-1,0 nohead lc 0" << endl
;
}
plot << "do for [loc_isid=0:"<<tilde_K.n_side()-1<<"] {" << endl
<< " loc_sid_ndof=" << degree()+1 << endl
<< " do for [i=0:loc_sid_ndof-1] {" << endl
<< " ip3=i+3" << endl
<< " splot gdat index loc_isid u 1:2:ip3 t sprintf(\"side %d: L%d\",loc_isid,i) w l lc 1" << endl
<< " pause pause_duration" << endl
<< " }" << endl
<< "}" << endl
<< "pause -1 \"<return>\"" << endl
;
plot.close();
// --------------------
// .gdat
// --------------------
ofstream gdat (gdat_name.c_str());
cerr << "! file \"" << gdat_name << "\" created" << endl;
filelist += " " + gdat_name;
gdat << setprecision(std::numeric_limits<T>::digits10)
<< "# basis " << base::name() << endl
<< "# element " << tilde_K.name() << endl
<< "# degree " << degree() << endl
;
Eigen::Matrix<T,Eigen::Dynamic,1> value;
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;
size_type loc_sid_ndof = first_sid_idof (tilde_K,loc_isid+1) - first_sid_idof (tilde_K,loc_isid);
gdat << "# side " << loc_isid << ": size " << loc_sid_ndof << endl
<< "# x y";
for (size_type loc_sid_idof = 0; loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
gdat << " L" << loc_sid_idof+1;
}
gdat << endl;
switch (hat_S.variant()) {
case reference_element::e: {
for (size_type i = 0; i <= nsub; i++) {
point_basic<T> sid_hat_x (T(int(i))/T(int(nsub)));
evaluate_on_side (tilde_K, sid, sid_hat_x, value);
point_basic<T> tilde_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
gdat << tilde_x[0] << " " << tilde_x[1];
for (size_type loc_idof = first_sid_idof (tilde_K,loc_isid), last_loc_idof = first_sid_idof (tilde_K,loc_isid+1);
loc_idof < last_loc_idof; ++loc_idof) {
gdat << " " << value[loc_idof];
}
gdat << endl;
}
gdat << endl << endl;
break;
}
default:
error_macro ("unexpected side type `"<<hat_S.name()<<"'");
}
}
// -----------
if (execute) {
// -----------
string command = "gnuplot " + plot_name;
if (verbose) clog << "! " << command << endl;
int status = system (command.c_str());
}
// -----------
if (clean) {
// -----------
string command = "/bin/rm -f " + filelist;
if (verbose) clog << "! " << command << endl;
int status = system (command.c_str());
}
#ifdef TODO
#endif // TODO
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class basis_fem_sides<T>;
_RHEOLEF_instanciation(Float)
}// namespace rheolef
basis_fem_sides.h
Float
see the Float page for the full documentation
rheolef::basis_basic
Definition: basis.h:532
rheolef::basis_fem_sides::basis_fem_sides
basis_fem_sides(const basis_basic< T > &side_basis)
Definition: basis_fem_sides.cc:71
rheolef::basis_fem_sides::_side_basis
basis_basic< T > _side_basis
Definition: basis_fem_sides.h:108
rheolef::basis_fem_sides::evaluate_on_side
void evaluate_on_side(reference_element hat_K, const side_information_type &sid, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
Definition: basis_fem_sides.cc:156
rheolef::basis_fem_sides::size_type
reference_element::size_type size_type
Definition: basis_fem_sides.h:54
rheolef::basis_fem_sides::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_sides.cc:188
rheolef::basis_fem_sides::degree
size_type degree() const
Definition: basis_fem_sides.h:66
rheolef::basis_fem_sides::_initialize_cstor_sizes
void _initialize_cstor_sizes() const
Definition: basis_fem_sides.cc:92
rheolef::basis_fem_sides::~basis_fem_sides
~basis_fem_sides()
Definition: basis_fem_sides.cc:67
rheolef::basis_fem_sides::local_ndof_on_side
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) const
Definition: basis_fem_sides.cc:178
rheolef::basis_fem_sides::put_scalar_valued
virtual void put_scalar_valued(std::ostream &os, reference_element hat_K) const
Definition: basis_fem_sides.cc:219
rheolef::basis_fem_sides::_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_sides.cc:206
rheolef::basis_fem_sides::family_index
size_type family_index() const
Definition: basis_fem_sides.h:65
rheolef::basis_fem_sides::hat_node
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
Definition: basis_fem_sides.cc:146
rheolef::basis_fem_sides::family_name
std::string family_name() const
Definition: basis_fem_sides.h:64
rheolef::basis_fem_sides::_initialize_data
void _initialize_data(reference_element hat_K) const
Definition: basis_fem_sides.cc:117
rheolef::basis_option::set_restricted_to_sides
void set_restricted_to_sides(bool r=true)
Definition: basis_option.h:281
rheolef::basis_rep
Definition: basis.h:209
rheolef::basis_rep::is_continuous
bool is_continuous() const
Definition: basis.h:240
rheolef::basis_rep::name
std::string name() const
Definition: basis.h:228
rheolef::basis_rep::size_type
reference_element::size_type size_type
Definition: basis.h:214
rheolef::basis_rep::_piola_fem
piola_fem< T > _piola_fem
Definition: basis.h:394
rheolef::basis_rep::_sopt
basis_option _sopt
Definition: basis.h:393
rheolef::basis_rep::option
const basis_option & option() const
Definition: basis.h:238
rheolef::basis_rep::standard_naming
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
rheolef::basis_rep::_name
std::string _name
Definition: basis.h:392
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::size
size_type size() const
Definition: reference_element.h:102
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::dimension
size_type dimension() const
Definition: reference_element.h:101
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::n_side
size_type n_side() const
Definition: reference_element.h:104
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::value
rheolef::std value
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
eigen_util.h
equispaced.icc
fekete.icc
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.clean
clean
Definition: mkgeo_ball.sh:335
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
mkgeo_ball.command
string command
Definition: mkgeo_ball.sh:136
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.verbose
verbose
Definition: mkgeo_ball.sh:133
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
mkgeo_couette.basename
string basename
Definition: mkgeo_couette.sh:73
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::evaluate_on_side
void evaluate_on_side(const geo_basic< float_type, M > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Matrix< Result, Eigen::Dynamic, 1 > &value) const
Definition: field_expr_recursive.h:894
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::reference_element_face_transformation
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
Definition: reference_element_face_transformation.cc:33
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
rheolef::side_information_type::loc_isid
size_t loc_isid
Definition: reference_element_face_transformation.h:38
warburton.icc