Loading...
Searching...
No Matches
field_expr_terminal.cc
Go to the documentation of this file.
//
//
// field_expr_terminal: terminals of nonn-linear expressions
// 1) class-function
// 1.1) base class for the class-function family
// 1.2) general function or class-function
// 1.3) normal to a surface
// 1.4) h_local
// 1.5) penalty
// 2) field and such
// 2.1) field
// 2.2) jump of a field
// 3) convected field, as compose(uh,X) where X is a characteristic
//
#include "rheolef/field_expr_terminal.h"
#include "rheolef/geo_domain.h"
#include "rheolef/piola_util.h"
#include "rheolef/field_evaluate.h"
#include "rheolef/memorized_disarray.h"
namespace rheolef { namespace details {
// ---------------------------------------------------------------------------
// 1.1) base class for all class-functions
// ---------------------------------------------------------------------------
template<class T>
field_expr_v2_nonlinear_terminal_function_base_rep<T>::field_expr_v2_nonlinear_terminal_function_base_rep ()
: _pops()
{
}
template<class T>
field_expr_v2_nonlinear_terminal_function_base_rep<T>::field_expr_v2_nonlinear_terminal_function_base_rep (const field_expr_v2_nonlinear_terminal_function_base_rep<T>& x)
: _pops(x._pops)
{
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_base_rep<T>::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
_pops = pops;
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_base_rep<T>::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
_pops = pops;
}
// TODO: move as geo::get_side (K, sid) ?
template<class T, class M>
const geo_element&
global_get_side (
const geo_basic<T,M>& omega_L,
const geo_element& L,
const side_information_type& sid)
{
typedef geo_element::size_type size_type;
size_type dis_isid = std::numeric_limits<size_type>::max();
switch (sid.dim) {
case 0: {
size_type dis_inod = L [sid.loc_isid];
dis_isid = omega_L.dis_inod2dis_iv (dis_inod);
break;
}
case 1: dis_isid = L.edge (sid.loc_isid); break;
case 2: dis_isid = L.face (sid.loc_isid); break;
default: error_macro ("domain: unexpected side dimension " << sid.dim);
}
// TODO: pas tres clair quand on est sur un geo_domain_indirect ou un geo_domain
// dans quel systeme de numerotation doit appartenir la face ?
if (omega_L.variant() != geo_abstract_base_rep<T>::geo_domain_indirect) {
return omega_L.dis_get_geo_element (sid.dim, dis_isid);
} else {
return omega_L.get_background_geo().dis_get_geo_element (sid.dim, dis_isid);
}
}
template<class T>
template<class M>
const geo_element&
field_expr_v2_nonlinear_terminal_function_base_rep<T>::get_side (
const geo_basic<T,M>& omega_L,
const geo_element& L,
const side_information_type& sid) const
{
return global_get_side (omega_L, L, sid);
}
// ----------------------------------------------------------------------------
// 1.3) normal() class-function
// ----------------------------------------------------------------------------
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const function_type&)
: base(),
_is_on_interface (false),
_is_inside_on_local_sides(false)
{
}
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >& x)
: base(x),
_is_on_interface (x._is_on_interface),
_is_inside_on_local_sides(x._is_inside_on_local_sides)
{
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (pops, iopt);
_is_on_interface = iopt._is_on_interface;
_is_inside_on_local_sides = iopt._is_inside_on_local_sides;
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (Xh, pops, iopt);
_is_on_interface = iopt._is_on_interface;
_is_inside_on_local_sides = iopt._is_inside_on_local_sides;
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::evaluate_internal (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
const T& sign,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = base::_pops.get_piola (omega_K, K);
size_type d = omega_K.dimension();
size_type loc_nnod = piola.size();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
const point_basic<float_type>& xi = piola[loc_inod].F;
value[loc_inod] = sign*normal_from_piola_transformation (omega_K, K, piola[loc_inod].DF, d);
}
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
evaluate_internal (omega_K, K, T(1), value);
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<normal_pseudo_function<T> >::evaluate_on_side (
const geo_basic<float_type,memory_type>& omega_L,
const geo_element& L,
const side_information_type& sid,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
if (_is_inside_on_local_sides) {
T sign = (L.dimension() > 1) ? sid.orient : ((sid.loc_isid == 0) ? -1 : 1);
evaluate_internal (omega_L, global_get_side(omega_L, L, sid), sign, value);
} else {
// BUG_INT_BDR_JUMP: check _is_on_interface and send sign of side in base::_omega
evaluate_internal (omega_L, global_get_side(omega_L, L, sid), T(1), value);
}
}
// ----------------------------------------------------------------------------
// 1.4) h_local() class-function
// ----------------------------------------------------------------------------
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const function_type&)
: base()
{
}
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >& x)
: base(x)
{
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (pops, iopt);
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (Xh, pops, iopt);
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = base::_pops.get_piola (omega_K, K);
size_type loc_nnod = piola.size();
value.resize (loc_nnod);
size_type K_map_d = K.dimension();
if (K_map_d > 0) {
// usual case: K is at least one edge
reference_element hat_K = K.variant();
T meas_hat_K = measure(hat_K);
size_type d = omega_K.dimension();
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
// TODO DVT_PIOLA_LOOP_ORDER_1: improve: when omega_K.order=1 then DF=constant and value[*]=constant
T det_DF_i = det_jacobian_piola_transformation (piola[loc_inod].DF, d, K_map_d);
value[loc_inod] = pow (abs(det_DF_i)/meas_hat_K, 1./K_map_d);
trace_macro("h_local("<<omega_K.name()<<",K="<<K.name()<<K.dis_ie()<<")="<< value[loc_inod]);
}
return;
}
// special case: d == 0 i.e. an "interface" between two edges in 1D
// then h_local = average of two edge length
size_type L_dis_ie[2] = {K.master(0), K.master(1)};
size_t n = (L_dis_ie[0] == std::numeric_limits<size_type>::max()) ? 0 :
(L_dis_ie[1] == std::numeric_limits<size_type>::max()) ? 1 : 2;
if (n == 0) {
// special 0d meshes:
// dis_nelt_1d=0: special 0d mesh, e.g. omega["sides"] when omega is 1d
// n=0: isolated point(s) in a 0D mesh
// => h_local=1 by default
value.fill (result_type(1));
trace_macro("h_local("<<omega_K.name()<<",K="<<K.name()<<K.dis_ie()<<")=1");
} else {
size_t dis_nelt_1d = omega_K.sizes().ownership_by_dimension[1].dis_size();
check_macro (dis_nelt_1d != 0, "h_local: invalid \""<<omega_K.name()<<"\" mesh with missing 1D connectivity");
size_t L_map_d = 1;
Float h_loc = 0;
for (size_type i = 0; i < n; ++i) {
const geo_element& Li = omega_K.dis_get_geo_element (L_map_d, L_dis_ie[i]);
const point_basic<float_type>& x0 = omega_K.dis_node (Li[0]);
const point_basic<float_type>& x1 = omega_K.dis_node (Li[1]);
h_loc += norm (x1-x0)/n;
}
value.fill (result_type(h_loc));
trace_macro("h_local("<<omega_K.name()<<",K="<<K.name()<<K.dis_ie()<<")="<< h_loc);
}
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >::evaluate_on_side (
const geo_basic<float_type,memory_type>& omega_L,
const geo_element& L,
const side_information_type& sid,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
evaluate (omega_L, global_get_side(omega_L, L, sid), value);
}
// ----------------------------------------------------------------------------
// 1.5) penalty() class-function
// ----------------------------------------------------------------------------
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const function_type&)
: base()
{
}
template<class T>
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::field_expr_v2_nonlinear_terminal_function_rep(
const field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >& x)
: base(x)
{
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (pops, iopt);
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
base::initialize (Xh, pops, iopt);
}
template<class T>
T
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::evaluate_measure (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K) const
{
const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = base::_pops.get_piola (omega_K, K);
size_type loc_nnod = piola.size();
size_type d = omega_K.dimension();
size_type map_d = K.dimension();
T meas_max_K = 0;
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
// TODO DVT_PIOLA_LOOP_ORDER_1: improve: when omega_K.order=1 then DF=constant and value[*]=constant
T meas_K = det_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
meas_max_K = std::max (meas_max_K, abs(meas_K));
}
trace_macro("meas("<<omega_K.name()<<",K="<<K.name()<<K.dis_ie()<<")="<< meas_max_K);
return meas_max_K;
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::evaluate_internal (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
const geo_element& L,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
// meas(partial L)/meas(L) ~= nside*meas(K)/meas(L) : TODO try to improve it for non-uniform meshes
size_type n_sides = L.n_subgeo (K.dimension());
T meas_side_L = evaluate_measure (omega_K, K);
T meas_L = evaluate_measure (omega_K, L);
T result = n_sides*meas_side_L/meas_L;
trace_macro("penalty("<<omega_K.name()<<",K="<<K.name()<<K.dis_ie()<<")="<< result);
const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = base::_pops.get_piola (omega_K, K);
size_type loc_nnod = piola.size();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] = result;
}
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
size_type L_map_d = K.dimension() + 1;
size_type L_dis_ie0 = K.master(0);
check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated side K="<<K.name()<<K.dis_ie()<< " in mesh \""<<omega_K.name()<<"\"");
const geo_element& L0 = omega_K.dis_get_geo_element (L_map_d, L_dis_ie0);
evaluate_internal (omega_K, K, L0, value);
size_type L_dis_ie1 = K.master(1);
if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
// K is a boundary side: nothing more to do
true;
} else {
// K is an internal side
const geo_element& L1 = omega_K.dis_get_geo_element (L_map_d, L_dis_ie1);
Eigen::Matrix<result_type,Eigen::Dynamic,1> value1;
evaluate_internal (omega_K, K, L0, value);
evaluate_internal (omega_K, K, L1, value1);
for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
value[loc_idof] = std::max (value [loc_idof], value1 [loc_idof]);
}
}
}
template<class T>
void
field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >::evaluate_on_side (
const geo_basic<float_type,memory_type>& omega_L,
const geo_element& L,
const side_information_type& sid,
Eigen::Matrix<result_type,Eigen::Dynamic,1>& value) const
{
evaluate (omega_L, global_get_side(omega_L, L, sid), value);
}
// ----------------------------------------------------------------------------
// 2) field and such
// ----------------------------------------------------------------------------
// 2.1) field
// ----------------------------------------------------------------------------
template<class T, class M, details::differentiate_option::type Diff>
void
field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::initialize (
const piola_on_pointset<T>& pops,
const integrate_option& iopt)
{
_uh.dis_dof_update();
_u_test.initialize (pops, iopt);
}
template<class T, class M, details::differentiate_option::type Diff>
void
field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::initialize (
const space_basic<T,M>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
check_macro (_uh.get_geo().get_background_geo().name() == Xh.get_geo().get_background_geo().name(),
"incompatible field on " << _uh.get_geo().get_background_geo().name() << " for evaluation on " << Xh.get_geo().get_background_geo().name());
_uh.dis_dof_update();
_u_test.initialize (Xh, pops, iopt);
}
template<class T, class M, details::differentiate_option::type Diff>
template<class Value>
void
field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<Value,Eigen::Dynamic,1>& value) const
{
Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> diff_phij_xi;
_u_test.template evaluate <Value,Diff> (omega_K, K, _gopt, diff_phij_xi);
field_evaluate_continued (_uh, omega_K, K, diff_phij_xi, value);
}
template<class T, class M, details::differentiate_option::type Diff>
template<class Value>
void
field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate_on_side (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
const side_information_type& sid,
Eigen::Matrix<Value,Eigen::Dynamic,1>& value) const
{
Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> diff_phij_xi;
bool do_local_component_assembly = true;
_u_test.template evaluate_on_side<Value,Diff> (omega_K, K, sid, _gopt, diff_phij_xi, do_local_component_assembly);
field_evaluate_continued (_uh, omega_K, K, diff_phij_xi, value);
}
template<class T, class M, details::differentiate_option::type Diff>
space_constant::valued_type
field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::valued_tag() const
{
switch (Diff) {
case details::differentiate_option::none: {
return _uh.valued_tag();
}
case details::differentiate_option::gradient: {
space_constant::valued_type uh_valued_tag = _uh.valued_tag();
switch (uh_valued_tag) {
case space_constant::scalar: return space_constant::vector;
case space_constant::vector: return space_constant::unsymmetric_tensor;
case space_constant::tensor:
case space_constant::unsymmetric_tensor:
return space_constant::tensor3;
default:
error_macro ("unsupported "<<uh_valued_tag<< "-valued field for the grad() operator");
return space_constant::last_valued;
}
}
case details::differentiate_option::divergence: {
space_constant::valued_type uh_valued_tag = _uh.valued_tag();
switch (uh_valued_tag) {
case space_constant::vector: return space_constant::scalar;
case space_constant::tensor:
case space_constant::unsymmetric_tensor:
return space_constant::vector;
default:
error_macro ("unsupported "<<uh_valued_tag<< "-valued field for the div() operator");
return space_constant::last_valued;
}
}
case details::differentiate_option::curl: {
space_constant::valued_type uh_valued_tag = _uh.valued_tag();
switch (uh_valued_tag) {
case space_constant::scalar: return space_constant::vector;
case space_constant::vector: {
size_type d = _uh.get_geo().dimension();
return (d==2) ? space_constant::scalar : space_constant::vector;
}
default:
error_macro ("unsupported "<<uh_valued_tag<< "-valued field for the curl() operator");
return space_constant::last_valued;
}
}
}
}
// ----------------------------------------------------------------------------
// 2.2) field jump
// ----------------------------------------------------------------------------
template<class T, class M>
void
field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
_expr0.initialize (pops, iopt);
_expr1.initialize (pops, iopt);
}
template<class T, class M>
void
field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
_expr0.initialize (Xh, pops, iopt);
_expr1.initialize (Xh, pops, iopt);
}
template<class T, class M>
template<class Result>
void
field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
{
size_type L_map_d = K.dimension() + 1;
size_type L_dis_ie0, L_dis_ie1;
side_information_type sid0, sid1;
L_dis_ie0 = K.master(0);
const geo_basic<T,M>& omega_L = omega_K.get_background_geo();
check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
const geo_element& L0 = omega_L.dis_get_geo_element (L_map_d, L_dis_ie0);
L0.get_side_informations (K, sid0);
_expr0.evaluate_on_side (omega_L, L0, sid0, value);
L_dis_ie1 = K.master(1);
if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
// K is a boundary side
// jump(v)=average(v)=inner(v)=outer(v)=v on the boundary
// TODO: DVT_EIGEN_BLAS2 : pb when Value=point
for (size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
for (size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
value(loc_inod,loc_jdof) = value(loc_inod,loc_jdof);
}}
} else {
// K is an internal side
Eigen::Matrix<Result,Eigen::Dynamic,1> value1;
const geo_element& L1 = omega_L.dis_get_geo_element (L_map_d, L_dis_ie1);
L1.get_side_informations (K, sid1);
_expr1.evaluate_on_side (omega_L, L1, sid1, value1);
for (size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
for (size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
value(loc_inod,loc_jdof) = _c0*value(loc_inod,loc_jdof) + _c1*value1(loc_inod,loc_jdof);
}}
}
}
// ---------------------------------------------------------------------------
// 3) contains a convected field, as compose(uh,X) where X is a characteristic
// ---------------------------------------------------------------------------
template<class T, class M, class Value>
void
interpolate_pass2_valued (
const field_basic<T,M>& uh,
const disarray<point_basic<T>,M>& x,
const disarray<index_set,M>& ie2dis_ix, // K -> list of ix
const disarray<point_basic<T>,M>& hat_y, // ix -> hat_y
disarray<Value,M>& ux);
template<class T, class M>
void
field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>::initialize (
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
// handle the case when omega (destination) is different from _X.get_displacement.get_geo (origin)
// as in mesh adaptation loops
_uh.dis_dof_update();
_fops.initialize (_uh.get_space().get_basis(), pops);
// coq is stored in _X : it is independent of _uh ; depend only of _uh.space and quadrature nodes (qopt)
const characteristic_on_quadrature<T,M> coq = _X.get_pre_computed (_uh.get_space(), _X.get_displacement(), pops);
const characteristic_on_quadrature_rep<T,M>& coq_r = coq.data();
switch (_uh.valued_tag()) {
case space_constant::scalar:
_scalar_val.resize (coq_r._yq.ownership());
interpolate_pass2_valued (_uh, coq_r._yq, coq_r._ie2dis_ix, coq_r._hat_y, _scalar_val);
break;
case space_constant::vector:
_vector_val.resize (coq_r._yq.ownership());
interpolate_pass2_valued (_uh, coq_r._yq, coq_r._ie2dis_ix, coq_r._hat_y, _vector_val);
break;
default: error_macro("unsupported "<<_uh.valued()<<"-valued charateristic");
}
_start_q = 0;
}
template<class T, class M>
void
field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>::initialize (
const space_basic<float_type,memory_type>& Xh,
const piola_on_pointset<float_type>& pops,
const integrate_option& iopt)
{
fatal_macro ("characteristic: not supported with interpolate (HINT: use integrate with Gauss-Lobatto quadrature)");
}
template<class T, class M>
template<class Value>
void
field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>::evaluate (
const geo_basic<float_type,memory_type>& omega_K,
const geo_element& K,
Eigen::Matrix<Value,Eigen::Dynamic,1>& value) const
{
const disarray<Value,M>& uq = memorized_disarray<T,Value>().get (*this);
typename disarray<Value>::const_iterator uq_K = uq.begin();
size_type loc_nnod = _fops.get_piola_on_pointset().get_basis_on_pointset().nnod(K);
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; loc_inod++, _start_q++) {
value [loc_inod] = uq_K [_start_q];
}
}
// ----------------------------------------------------------------------------
// instanciation in library (v2)
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation_base_seq_only(T) \
template class field_expr_v2_nonlinear_terminal_function_base_rep<T>; \
template class field_expr_v2_nonlinear_terminal_function_rep< normal_pseudo_function<T> >; \
template class field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >; \
template class field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >; \
// extra when Result=tensor3, tensor4
#define _RHEOLEF_instanciation_base_both_members_diff(T,M,Result,Diff) \
template void field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate ( \
const geo_basic<T,M>& omega_K, \
const geo_element& K, \
Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
template void field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate_on_side ( \
const geo_basic<T,M>& omega_L, \
const geo_element& L, \
const side_information_type& sid, \
Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
#define _RHEOLEF_instanciation_base_both_members_alls(T,M,Result) \
_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::none) \
_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::gradient) \
_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::divergence) \
_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::curl) \
template void field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>::evaluate ( \
const geo_basic<float_type,memory_type>& omega_K, \
const geo_element& K, \
Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
template void field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>::evaluate ( \
const geo_basic<float_type,memory_type>& omega_K, \
const geo_element& K, \
Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
#define _RHEOLEF_instanciation_base_both(T,M) \
template const geo_element& global_get_side ( \
const geo_basic<T,M>&, const geo_element&, const side_information_type&); \
template const geo_element& field_expr_v2_nonlinear_terminal_function_base_rep<T>::get_side ( \
const geo_basic<T,M>& omega_L, \
const geo_element& L, \
const side_information_type& sid) const; \
template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::none>; \
template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::gradient>; \
template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::divergence>; \
template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::curl>; \
template class field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>; \
template class field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>; \
_RHEOLEF_instanciation_base_both_members_alls(T,M,T) \
_RHEOLEF_instanciation_base_both_members_alls(T,M,point_basic<T>) \
_RHEOLEF_instanciation_base_both_members_alls(T,M,tensor_basic<T>) \
_RHEOLEF_instanciation_base_both_members_alls(T,M,tensor3_basic<T>) \
_RHEOLEF_instanciation_base_both_members_alls(T,M,tensor4_basic<T>) \
#define _RHEOLEF_instanciation_seq(T) \
_RHEOLEF_instanciation_base_seq_only(T) \
_RHEOLEF_instanciation_base_both(T,sequential) \
#define _RHEOLEF_instanciation_dis(T) \
_RHEOLEF_instanciation_base_both(T,distributed) \
_RHEOLEF_instanciation_seq(Float)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation_dis(Float)
#endif // _RHEOLEF_HAVE_MPI
}} // namespace rheolef::details
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
rheolef::characteristic_on_quadrature
Definition characteristic.h:150
rheolef::details::field_expr_v2_nonlinear_terminal_field_dg_rep::initialize
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
Definition field_expr_terminal.cc:500
rheolef::details::field_expr_v2_nonlinear_terminal_field_dg_rep::evaluate
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
rheolef::details::field_expr_v2_nonlinear_terminal_field_dg_rep::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:1127
rheolef::details::field_expr_v2_nonlinear_terminal_field_o_characteristic_rep::initialize
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
Definition field_expr_terminal.cc:570
rheolef::details::field_expr_v2_nonlinear_terminal_field_o_characteristic_rep::evaluate
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition field_expr_terminal.cc:607
rheolef::details::field_expr_v2_nonlinear_terminal_field_o_characteristic_rep::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:1269
rheolef::details::field_expr_v2_nonlinear_terminal_field_rep::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:636
rheolef::details::field_expr_v2_nonlinear_terminal_field_rep::initialize
void initialize(const piola_on_pointset< T > &pops, const integrate_option &iopt)
Definition field_expr_terminal.cc:402
rheolef::details::field_expr_v2_nonlinear_terminal_field_rep::evaluate
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition field_expr_terminal.cc:424
rheolef::details::field_expr_v2_nonlinear_terminal_field_rep::evaluate_on_side
void evaluate_on_side(const geo_basic< float_type, memory_type > &omega_L, const geo_element &L, const side_information_type &sid, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition field_expr_terminal.cc:436
rheolef::details::field_expr_v2_nonlinear_terminal_field_rep::valued_tag
space_constant::valued_type valued_tag() const
Definition field_expr_terminal.cc:449
rheolef::details::field_expr_v2_nonlinear_terminal_function_base_rep
Definition field_expr_terminal.h:101
rheolef::details::field_expr_v2_nonlinear_terminal_function_base_rep::initialize
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
Definition field_expr_terminal.cc:57
rheolef::details::field_expr_v2_nonlinear_terminal_function_base_rep::get_side
const geo_element & get_side(const geo_basic< float_type, M > &omega_K, const geo_element &K, const side_information_type &sid) const
Definition field_expr_terminal.cc:103
rheolef::details::field_expr_v2_nonlinear_terminal_function_base_rep::field_expr_v2_nonlinear_terminal_function_base_rep
field_expr_v2_nonlinear_terminal_function_base_rep()
Definition field_expr_terminal.cc:46
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< h_local_pseudo_function< T > >::base
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
Definition field_expr_terminal.h:456
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< h_local_pseudo_function< T > >::result_type
T result_type
Definition field_expr_terminal.h:460
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< h_local_pseudo_function< T > >::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:457
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< normal_pseudo_function< T > >::base
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
Definition field_expr_terminal.h:360
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< normal_pseudo_function< T > >::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:361
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< penalty_pseudo_function< T > >::base
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
Definition field_expr_terminal.h:544
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep< penalty_pseudo_function< T > >::size_type
geo_element::size_type size_type
Definition field_expr_terminal.h:545
rheolef::details::field_expr_v2_nonlinear_terminal_function_rep
Definition field_expr_terminal.h:144
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::disarray::const_iterator
rep::base::const_iterator const_iterator
Definition disarray.h:503
rheolef::field_basic
Definition field.h:219
rheolef::geo_abstract_base_rep
abstract base interface class
Definition geo.h:248
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_element::size_type
reference_element::size_type size_type
Definition geo_element.h:125
rheolef::geo_element::dimension
size_type dimension() const
Definition geo_element.h:167
rheolef::geo_element::master
size_type master(bool i) const
Definition geo_element.h:165
rheolef::geo_element::variant
variant_type variant() const
Definition geo_element.h:161
rheolef::geo_element::get_side_informations
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
Definition geo_element.cc:185
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition geo_element.h:163
rheolef::geo_element::name
char name() const
Definition geo_element.h:169
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition integrate_option.h:125
rheolef::integrate_option::_is_on_interface
bool _is_on_interface
Definition integrate_option.h:180
rheolef::integrate_option::_is_inside_on_local_sides
bool _is_inside_on_local_sides
Definition integrate_option.h:180
rheolef::piola_on_pointset
Definition piola_on_pointset.h:142
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::smart_pointer_base::data
const T & data() const
Definition smart_pointer.h:266
rheolef::space_basic
the finite element space
Definition space.h:382
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
measure
const Float measure
Definition edge.icc:65
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_instanciation_seq
#define _RHEOLEF_instanciation_seq(T)
Definition field_expr_terminal.cc:674
_RHEOLEF_instanciation_dis
#define _RHEOLEF_instanciation_dis(T)
Definition field_expr_terminal.cc:678
rheolef::details::global_get_side
const geo_element & global_get_side(const geo_basic< T, M > &omega_L, const geo_element &L, const side_information_type &sid)
Definition field_expr_terminal.cc:75
rheolef::details::interpolate_pass2_valued
void interpolate_pass2_valued(const field_basic< T, M > &uh, const disarray< point_basic< T >, M > &x, const disarray< index_set, M > &ie2dis_ix, const disarray< point_basic< T >, M > &hat_y, disarray< Value, M > &ux)
Definition interpolate.cc:122
rheolef::space_constant::valued_type
valued_type
Definition space_constant.h:133
rheolef::space_constant::tensor
@ tensor
Definition space_constant.h:136
rheolef::space_constant::scalar
@ scalar
Definition space_constant.h:134
rheolef::space_constant::unsymmetric_tensor
@ unsymmetric_tensor
Definition space_constant.h:137
rheolef::space_constant::vector
@ vector
Definition space_constant.h:135
rheolef::space_constant::tensor3
@ tensor3
Definition space_constant.h:138
rheolef::space_constant::last_valued
@ last_valued
Definition space_constant.h:141
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::normal_from_piola_transformation
point_basic< T > normal_from_piola_transformation(const geo_basic< T, M > &omega, const geo_element &S, const tensor_basic< T > &DF, size_t d)
Definition piola_util.cc:229
rheolef::initialize
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
Definition field_expr_recursive.h:627
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
rheolef::field_evaluate_continued
void field_evaluate_continued(const field_basic< T, M > &uh, const geo_basic< T, M > &omega_K, const geo_element &K, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &phij_xi, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value)
Definition field_evaluate.cc:217
rheolef::evaluate
void evaluate(const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, 1 > &value) const
Definition field_expr_recursive.h:915
rheolef::det_jacobian_piola_transformation
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
Definition piola_util.cc:137
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:936
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
d
size_t d
Definition neumann-laplace-lambda.cc:28
rheolef::characteristic_on_quadrature_rep
Definition characteristic.h:118
rheolef::characteristic_on_quadrature_rep::_ie2dis_ix
disarray< index_set, M > _ie2dis_ix
Definition characteristic.h:145
rheolef::characteristic_on_quadrature_rep::_hat_y
disarray< point_basic< T >, M > _hat_y
Definition characteristic.h:146
rheolef::characteristic_on_quadrature_rep::_yq
disarray< point_basic< T >, M > _yq
Definition characteristic.h:147
rheolef::details::differentiate_option::divergence
@ divergence
Definition piola.h:45
rheolef::details::differentiate_option::gradient
@ gradient
Definition piola.h:44
rheolef::details::differentiate_option::none
@ none
Definition piola.h:43
rheolef::details::differentiate_option::curl
@ curl
Definition piola.h:46
rheolef::details::h_local_pseudo_function
Definition field_expr_terminal.h:443
rheolef::details::memorized_disarray
Definition memorized_disarray.h:44
rheolef::details::normal_pseudo_function
Definition field_expr_terminal.h:348
rheolef::details::penalty_pseudo_function
Definition field_expr_terminal.h:531
rheolef::piola
Definition piola.h:67
rheolef::piola::F
point_basic< T > F
Definition piola.h:79
rheolef::side_information_type
Definition reference_element_face_transformation.h:37
rheolef::side_information_type::dim
size_t dim
Definition reference_element_face_transformation.h:41
rheolef::side_information_type::orient
geo_element_indirect::orientation_type orient
Definition reference_element_face_transformation.h:40
rheolef::side_information_type::loc_isid
size_t loc_isid
Definition reference_element_face_transformation.h:38
M
Expr1::memory_type M
Definition vec_expr_v2.h:416