34#include "rheolef/field_expr_terminal.h"
35#include "rheolef/geo_domain.h"
36#include "rheolef/piola_util.h"
37#include "rheolef/field_evaluate.h"
38#include "rheolef/memorized_disarray.h"
40namespace rheolef {
namespace details {
73template<
class T,
class M>
81 size_type dis_isid = std::numeric_limits<size_type>::max();
85 dis_isid = omega_L.dis_inod2dis_iv (dis_inod);
88 case 1: dis_isid = L.edge (sid.
loc_isid);
break;
89 case 2: dis_isid = L.face (sid.
loc_isid);
break;
90 default:
error_macro (
"domain: unexpected side dimension " << sid.
dim);
95 return omega_L.dis_get_geo_element (sid.
dim, dis_isid);
97 return omega_L.get_background_geo().dis_get_geo_element (sid.
dim, dis_isid);
117 _is_on_interface (false),
118 _is_inside_on_local_sides(false)
125 _is_on_interface (x._is_on_interface),
126 _is_inside_on_local_sides(x._is_inside_on_local_sides)
135 base::initialize (pops, iopt);
146 base::initialize (Xh, pops, iopt);
156 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
158 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>&
piola = base::_pops.get_piola (omega_K, K);
161 value.resize (loc_nnod);
162 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
172 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
174 evaluate_internal (omega_K, K,
T(1), value);
182 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
184 if (_is_inside_on_local_sides) {
185 T sign = (L.dimension() > 1) ? sid.
orient : ((sid.
loc_isid == 0) ? -1 : 1);
186 evaluate_internal (omega_L,
global_get_side(omega_L, L, sid), sign, value);
213 base::initialize (pops, iopt);
222 base::initialize (Xh, pops, iopt);
229 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
231 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>&
piola = base::_pops.get_piola (omega_K, K);
233 value.resize (loc_nnod);
240 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
243 value[loc_inod] =
pow (abs(det_DF_i)/meas_hat_K, 1./K_map_d);
251 size_t n = (L_dis_ie[0] == std::numeric_limits<size_type>::max()) ? 0 :
252 (L_dis_ie[1] == std::numeric_limits<size_type>::max()) ? 1 : 2;
261 size_t dis_nelt_1d = omega_K.sizes().ownership_by_dimension[1].dis_size();
262 check_macro (dis_nelt_1d != 0,
"h_local: invalid \""<<omega_K.name()<<
"\" mesh with missing 1D connectivity");
266 const geo_element& Li = omega_K.dis_get_geo_element (L_map_d, L_dis_ie[i]);
269 h_loc +=
norm (x1-x0)/n;
281 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
306 base::initialize (pops, iopt);
315 base::initialize (Xh, pops, iopt);
323 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>&
piola = base::_pops.get_piola (omega_K, K);
328 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
331 meas_max_K = std::max (meas_max_K, abs(meas_K));
342 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
346 T meas_side_L = evaluate_measure (omega_K, K);
347 T meas_L = evaluate_measure (omega_K, L);
348 T result = n_sides*meas_side_L/meas_L;
350 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>&
piola = base::_pops.get_piola (omega_K, K);
352 value.resize (loc_nnod);
353 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
354 value[loc_inod] = result;
362 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
366 check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
367 "unexpected isolated side K="<<K.
name()<<K.
dis_ie()<<
" in mesh \""<<omega_K.name()<<
"\"");
368 const geo_element& L0 = omega_K.dis_get_geo_element (L_map_d, L_dis_ie0);
369 evaluate_internal (omega_K, K, L0, value);
371 if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
376 const geo_element& L1 = omega_K.dis_get_geo_element (L_map_d, L_dis_ie1);
377 Eigen::Matrix<result_type,Eigen::Dynamic,1> value1;
378 evaluate_internal (omega_K, K, L0, value);
379 evaluate_internal (omega_K, K, L1, value1);
380 for (
size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
381 value[loc_idof] = std::max (value [loc_idof], value1 [loc_idof]);
391 Eigen::Matrix<result_type,Eigen::Dynamic,1>& value)
const
400template<
class T,
class M, details::differentiate_option::type Diff>
406 _uh.dis_dof_update();
407 _u_test.initialize (pops, iopt);
409template<
class T,
class M, details::differentiate_option::type Diff>
416 check_macro (_uh.get_geo().get_background_geo().name() == Xh.get_geo().get_background_geo().name(),
417 "incompatible field on " << _uh.get_geo().get_background_geo().name() <<
" for evaluation on " << Xh.get_geo().get_background_geo().name());
418 _uh.dis_dof_update();
419 _u_test.initialize (Xh, pops, iopt);
421template<
class T,
class M, details::differentiate_option::type Diff>
427 Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
const
429 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> diff_phij_xi;
430 _u_test.template evaluate <Value,Diff> (omega_K, K, _gopt, diff_phij_xi);
433template<
class T,
class M, details::differentiate_option::type Diff>
440 Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
const
442 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> diff_phij_xi;
443 bool do_local_component_assembly =
true;
444 _u_test.template evaluate_on_side<Value,Diff> (omega_K, K, sid, _gopt, diff_phij_xi, do_local_component_assembly);
447template<
class T,
class M, details::differentiate_option::type Diff>
453 return _uh.valued_tag();
457 switch (uh_valued_tag) {
464 error_macro (
"unsupported "<<uh_valued_tag<<
"-valued field for the grad() operator");
470 switch (uh_valued_tag) {
476 error_macro (
"unsupported "<<uh_valued_tag<<
"-valued field for the div() operator");
482 switch (uh_valued_tag) {
489 error_macro (
"unsupported "<<uh_valued_tag<<
"-valued field for the curl() operator");
498template<
class T,
class M>
504 _expr0.initialize (pops, iopt);
505 _expr1.initialize (pops, iopt);
507template<
class T,
class M>
514 _expr0.initialize (Xh, pops, iopt);
515 _expr1.initialize (Xh, pops, iopt);
517template<
class T,
class M>
518template<
class Result>
523 Eigen::Matrix<Result,Eigen::Dynamic,1>& value)
const
530 check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
531 "unexpected isolated mesh side K.dis_ie="<<K.
dis_ie());
532 const geo_element& L0 = omega_L.dis_get_geo_element (L_map_d, L_dis_ie0);
534 _expr0.evaluate_on_side (omega_L, L0, sid0, value);
536 if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
540 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
541 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
542 value(loc_inod,loc_jdof) = value(loc_inod,loc_jdof);
546 Eigen::Matrix<Result,Eigen::Dynamic,1> value1;
547 const geo_element& L1 = omega_L.dis_get_geo_element (L_map_d, L_dis_ie1);
549 _expr1.evaluate_on_side (omega_L, L1, sid1, value1);
550 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
551 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
552 value(loc_inod,loc_jdof) = _c0*value(loc_inod,loc_jdof) + _c1*value1(loc_inod,loc_jdof);
559template<
class T,
class M,
class Value>
568template<
class T,
class M>
576 _uh.dis_dof_update();
577 _fops.initialize (_uh.get_space().get_basis(), pops);
582 switch (_uh.valued_tag()) {
584 _scalar_val.resize (coq_r.
_yq.ownership());
588 _vector_val.resize (coq_r.
_yq.ownership());
591 default:
error_macro(
"unsupported "<<_uh.valued()<<
"-valued charateristic");
595template<
class T,
class M>
602 fatal_macro (
"characteristic: not supported with interpolate (HINT: use integrate with Gauss-Lobatto quadrature)");
604template<
class T,
class M>
610 Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
const
614 size_type loc_nnod = _fops.get_piola_on_pointset().get_basis_on_pointset().nnod(K);
615 value.resize (loc_nnod);
616 for (
size_type loc_inod = 0; loc_inod < loc_nnod; loc_inod++, _start_q++) {
617 value [loc_inod] = uq_K [_start_q];
623#define _RHEOLEF_instanciation_base_seq_only(T) \
624template class field_expr_v2_nonlinear_terminal_function_base_rep<T>; \
625template class field_expr_v2_nonlinear_terminal_function_rep< normal_pseudo_function<T> >; \
626template class field_expr_v2_nonlinear_terminal_function_rep<h_local_pseudo_function<T> >; \
627template class field_expr_v2_nonlinear_terminal_function_rep<penalty_pseudo_function<T> >; \
630#define _RHEOLEF_instanciation_base_both_members_diff(T,M,Result,Diff) \
631template void field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate ( \
632 const geo_basic<T,M>& omega_K, \
633 const geo_element& K, \
634 Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
635template void field_expr_v2_nonlinear_terminal_field_rep<T,M,Diff>::evaluate_on_side ( \
636 const geo_basic<T,M>& omega_L, \
637 const geo_element& L, \
638 const side_information_type& sid, \
639 Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
641#define _RHEOLEF_instanciation_base_both_members_alls(T,M,Result) \
642_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::none) \
643_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::gradient) \
644_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::divergence) \
645_RHEOLEF_instanciation_base_both_members_diff(T,M,Result,details::differentiate_option::curl) \
646template void field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>::evaluate ( \
647 const geo_basic<float_type,memory_type>& omega_K, \
648 const geo_element& K, \
649 Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
650template void field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>::evaluate ( \
651 const geo_basic<float_type,memory_type>& omega_K, \
652 const geo_element& K, \
653 Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const; \
655#define _RHEOLEF_instanciation_base_both(T,M) \
656template const geo_element& global_get_side ( \
657 const geo_basic<T,M>&, const geo_element&, const side_information_type&); \
658template const geo_element& field_expr_v2_nonlinear_terminal_function_base_rep<T>::get_side ( \
659 const geo_basic<T,M>& omega_L, \
660 const geo_element& L, \
661 const side_information_type& sid) const; \
662template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::none>; \
663template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::gradient>; \
664template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::divergence>; \
665template class field_expr_v2_nonlinear_terminal_field_rep<T,M,details::differentiate_option::curl>; \
666template class field_expr_v2_nonlinear_terminal_field_dg_rep<T,M>; \
667template class field_expr_v2_nonlinear_terminal_field_o_characteristic_rep<T,M>; \
668 _RHEOLEF_instanciation_base_both_members_alls(T,M,T) \
669 _RHEOLEF_instanciation_base_both_members_alls(T,M,point_basic<T>) \
670 _RHEOLEF_instanciation_base_both_members_alls(T,M,tensor_basic<T>) \
671 _RHEOLEF_instanciation_base_both_members_alls(T,M,tensor3_basic<T>) \
672 _RHEOLEF_instanciation_base_both_members_alls(T,M,tensor4_basic<T>) \
674#define _RHEOLEF_instanciation_seq(T) \
675 _RHEOLEF_instanciation_base_seq_only(T) \
676 _RHEOLEF_instanciation_base_both(T,sequential) \
678#define _RHEOLEF_instanciation_dis(T) \
679 _RHEOLEF_instanciation_base_both(T,distributed) \
682#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
geo_element::size_type size_type
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
geo_element::size_type size_type
geo_element::size_type size_type
void initialize(const piola_on_pointset< T > &pops, const integrate_option &iopt)
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
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
space_constant::valued_type valued_tag() const
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
const geo_element & get_side(const geo_basic< float_type, M > &omega_K, const geo_element &K, const side_information_type &sid) const
field_expr_v2_nonlinear_terminal_function_base_rep()
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
geo_element::size_type size_type
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
geo_element::size_type size_type
field_expr_v2_nonlinear_terminal_function_base_rep< T > base
geo_element::size_type size_type
see the disarray page for the full documentation
rep::base::const_iterator const_iterator
abstract base interface class
generic mesh with rerefence counting
see the geo_element page for the full documentation
reference_element::size_type size_type
size_type dimension() const
size_type master(bool i) const
variant_type variant() const
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
see the integrate_option page for the full documentation
bool _is_inside_on_local_sides
see the reference_element page for the full documentation
#define trace_macro(message)
#define error_macro(message)
#define fatal_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_instanciation_seq(T)
#define _RHEOLEF_instanciation_dis(T)
const geo_element & global_get_side(const geo_basic< T, M > &omega_L, const geo_element &L, const side_information_type &sid)
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)
This file is part of Rheolef.
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)
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
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)
void evaluate(const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, 1 > &value) const
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
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
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
disarray< index_set, M > _ie2dis_ix
disarray< point_basic< T >, M > _hat_y
disarray< point_basic< T >, M > _yq