1# ifndef _RHEOLEF_FIELD_LAZY_TERMINAL_H
2# define _RHEOLEF_FIELD_LAZY_TERMINAL_H
37#include "rheolef/field_lazy.h"
52template<
class Derived>
53class field_lazy_base {
59 Derived&
derived() {
return *
static_cast< Derived*
>(
this); }
60 const Derived&
derived()
const {
return *
static_cast<const Derived*
>(
this); }
76template<
class T,
class M>
114template<
class T,
class M>
118template<
class T,
class M>
122 _uh.dis_dof_update();
124template<
class T,
class M>
131 std::vector<size_type> dis_idof;
132 _uh.get_space().get_constitution().assembly_dis_idof (omega_K, K, dis_idof);
134 uk.resize (loc_ndof);
135 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
137 uk [loc_jdof] = _uh.dis_dof (dis_jdof);
258 return get_space().get_constitution().get_geo();
265 _iopt._is_on_interface =
false;
266 _iopt._is_inside_on_local_sides =
false;
267 _expr.initialize (_domain, _iopt);
268 _prev_omega_K = omega_K;
269 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
278 _is_on_domain.resize (omega_K.geo_element_ownership(K_map_d));
279 std::fill (_is_on_domain.begin(), _is_on_domain.end(),
false);
280 if (_domain.map_dimension() == omega_K.map_dimension()) {
285 trace_macro (
"lazy_int(init): domain="<<_domain.name()<<
": same numbering as "<<omega_K.name()<<
"...");
286 size_type first_dis_ie = omega_K.geo_element_ownership (omega_K.map_dimension()).first_index();
287 trace_macro (
"lazy_int(init): first_dis_ie="<<first_dis_ie);
288 trace_macro (
"lazy_int(init): _is_on_domain.size="<<_is_on_domain.size());
289 for (
size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
294 check_macro (first_dis_ie <= dis_ie,
"invalid index");
297 check_macro (ie <= _is_on_domain.size(),
"invalid index");
298 _is_on_domain [ie] =
true;
301 trace_macro (
"lazy_int(init): domain="<<_domain.name()<<
": same numbering as "<<omega_K.name()<<
" done");
303 error_macro (
"geo_domain="<<_domain.name()<<
" subset of "<<omega_K.name()<<
": not yet");
305 }
else if (_domain.map_dimension()+1 == omega_K.map_dimension()) {
306 check_macro (_domain.get_background_geo() == omega_K,
"unexpected domain \""
307 <<_domain.name()<<
"\" as subset of \""<<omega_K.name()<<
"\"");
309 trace_macro (
"lazy_int(init): domain(bdry or sides)="<<_domain.name()<<
" subset of "<<omega_K.name()<<
"...");
310 for (
size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
316 check_macro (K_dis_ie != std::numeric_limits<size_type>::max(),
317 "unexpected isolated side S="<<S.
name()<<S.
dis_ie());
318 const geo_element& K = omega_K.dis_get_geo_element (K_map_d, K_dis_ie);
320 _is_on_domain.dis_entry (dis_ie) =
true;
322 _is_on_domain.dis_entry_assembly();
323 trace_macro (
"lazy_int(init): domain(bdry or sides)="<<_domain.name()<<
" subset of "<<omega_K.name()<<
" done");
325 error_macro (
"geo_domain(bdry or sides)="<<_domain.name()<<
" subset of "<<omega_K.name()<<
": not yet");
328 error_macro (
"unsupported domain \"" << _domain.name()
329 <<
"\" with map_dimension=" << _domain.map_dimension()
330 <<
" when integrated in geometry \"" << omega_K.name()
331 <<
"\" with map_dimension=" << omega_K.map_dimension());
369 std::array<std::vector<size_t>,
379 trace_macro (
"compute_idof_S2idof_K: sid_dim="<<sid_dim);
380 ndof_K = b.ndof (hat_K);
386 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
387 is_in_S[loc_isid].resize (ndof_K);
388 std::fill (is_in_S[loc_isid].begin(), is_in_S[loc_isid].end(),
false);
390 trace_macro (
"compute_idof_S2idof_K: loc_isid="<<loc_isid<<
", hat_S="<<hat_S.
name());
392 for (
size_type sub_sid_dim = 0; sub_sid_dim <= sid_dim; ++sub_sid_dim) {
393 size_type first_idof_K = b.first_idof_by_dimension (hat_K, sub_sid_dim);
394 trace_macro (
"compute_idof_S2idof_K: first_idof_K(hat_K="<<hat_K.
name()<<
",sub_sid_dim="<<sub_sid_dim<<
") = " <<first_idof_K);
395 for (
size_type loc_isub_sid = 0; loc_isub_sid < hat_S.
n_subgeo(sub_sid_dim); ++loc_isub_sid) {
399 switch (sub_sid_dim) {
406 loc_T_idx = loc_isid;
408 error_macro(
"compute_idof_S2idof_K: dof on edge in 3D: not yet");
415 loc_T_idx = loc_isid;
417 error_macro(
"compute_idof_S2idof_K: dof on face in 3D: not yet");
422 trace_macro (
"compute_idof_S2idof_K: ndof_K_on_T(hat_K.dim="<<hat_K.
dimension()<<
",hat_T="<<hat_T.
name()<<
") = " <<ndof_K_on_T);
423 trace_macro (
"compute_idof_S2idof_K: T_idx="<<loc_T_idx);
424 size_type start_idof_K = first_idof_K + ndof_K_on_T*loc_T_idx;
425 size_type last_idof_K = first_idof_K + ndof_K_on_T*(loc_T_idx+1);
426 for (
size_type idof_K = start_idof_K; idof_K < last_idof_K; ++idof_K) {
427 if (is_in_S [loc_isid] [idof_K])
continue;
428 is_in_S [loc_isid] [idof_K] =
true;
429 idof_S2idof_K_map [loc_isid] [idof_S] = idof_K;
430 trace_macro (
"idof_S2idof_K_map(isid="<<loc_isid<<
",idof_S="<<idof_S<<
") = "<<idof_K);
440 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
441 ndof_S [loc_isid] = 0;
442 for (
size_type idof_K = 0; idof_K < ndof_K; ++idof_K) {
443 if (is_in_S [loc_isid] [idof_K]) {
447 trace_macro (
"compute_idof_S2idof_K: ndof_S[isid="<<loc_isid<<
"]="<<ndof_S[loc_isid]);
453 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
454 ndof_S [loc_isid] = idof_S2idof_K_map [loc_isid].size();
455 idof_S2idof_K[loc_isid].resize (ndof_S [loc_isid]);
456 for (
auto p : idof_S2idof_K_map [loc_isid]) {
459 if (is_in_S [loc_isid] [idof_K]) {
460 idof_S2idof_K [loc_isid] [idof_S] = idof_K;
461 trace_macro (
"idof_K(isid="<<loc_isid<<
",idof_S="<<idof_S<<
") = "<< idof_K);
475 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.
dis_ie()) {
479 if (_domain.map_dimension() == omega_K.map_dimension()) {
480 trace_macro (
"lazy_int(eval): domain="<<_domain.name()<<
": same numbering as "<<omega_K.name()<<
"...");
481 _expr.evaluate (omega_K, K, uk);
482 trace_macro (
"lazy_int(eval): domain="<<_domain.name()<<
": same numbering as "<<omega_K.name()<<
" done");
484 trace_macro (
"lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<
" subset of "<<omega_K.name()<<
"...");
486 std::array<size_type, reference_element::max_side_by_variant> ndof_S;
488 check_macro (get_space().size() == 0,
"lazy_int(eval): "<<get_space().size()<<
"-component space: not yet supported");
490 uk = vector_element_type::Zero (ndof_K, 1);
491 trace_macro (
"lazy_int(eval): uk.size="<<uk.size());
495 for (
size_type loc_isid = 0, nsid = K.
n_subgeo(sid_dim); loc_isid < nsid; ++loc_isid) {
497 const geo_element& S = omega_K.dis_get_geo_element (sid_dim, dis_isid);
500 _expr.evaluate (_domain, S, us);
501 trace_macro (
"lazy_int(eval): us(loc_isid="<<loc_isid<<
").size="<<us.size());
503 for (
size_type idof_S = 0; idof_S < ndof_S [loc_isid]; ++idof_S) {
504 size_type idof_K = idof_S2idof_K[loc_isid][idof_S];
505 trace_macro (
"lazy_int(eval): uk[idof_K="<<idof_K<<
"] += us[idof_S="<<idof_S<<
"] = " << us[idof_S]);
506 uk[idof_K] += us[idof_S];
509 trace_macro (
"lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<
" subset of "<<omega_K.name()<<
" done");
512 _prev_omega_K = omega_K;
513 _prev_K_dis_ie = K.
dis_ie();
540 :
base1(new_macro(
rep(domain,expr,iopt))),
555 :
base1(new_macro(
rep(domname,expr,iopt))),
563 bool is_on_band()
const {
return base1::data().is_on_band(); }
572 { base1::data().evaluate (omega_K, K, uk); }
585 ,details::field_lazy_terminal_integrate <Expr>
589 const typename Expr::geo_type& domain,
599 details::is_field_expr_v2_variational_arg<Expr>::value
600 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
616 details::is_field_expr_quadrature_arg<Expr>::value
617 ,details::field_lazy_terminal_integrate <Expr>
630 details::is_field_expr_v2_variational_arg<Expr>::value
631 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
646 details::is_field_expr_quadrature_arg<Expr>::value
647 ,details::field_lazy_terminal_integrate <Expr>
651 const std::string& domname,
661 details::is_field_expr_v2_variational_arg<Expr>::value
662 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
666 const std::string& domname,
740 const geo_type& bgd_omega = get_space().get_constitution().get_background_geo();
742 "integrate on band: unexpected integration domain");
744 "do_integrate: incompatible integration domain "<<_gh.level_set().name() <<
" and test function based domain "
745 << bgd_omega.name());
746 _expr.initialize (_gh, _iopt);
747 _prev_omega_K = omega_K;
748 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
757 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.
dis_ie()) {
761 _expr.evaluate (omega_K, K, uk);
763 _prev_omega_K = omega_K;
764 _prev_K_dis_ie = K.
dis_ie();
797 bool is_on_band()
const {
return base1::data().is_on_band(); }
806 { base1::data().evaluate (omega_K, K, uk); }
817 ,details::field_lazy_terminal_integrate_band <Expr>
822 typename Expr::memory_type>&
gh,
831 details::is_field_expr_v2_variational_arg<Expr>::value
832 ,details::field_lazy_terminal_integrate_band <details::field_expr_quadrature_on_element<Expr>>
837 typename Expr::memory_type>&
gh,
860 using geo_type = geo_basic <float_type,memory_type>;
887 template <
class Value>
888 typename std::enable_if<
897 template <
class Value>
898 typename std::enable_if<
935 "interpolate: incompatible space \""<<_Xh.name()<<
"\" with geometry \""
938 const piola_fem<float_type>& pf = b.get_piola_fem();
940 _pops.initialize (omega_K.get_piola_basis(), b, iopt);
941 _expr.initialize (_Xh, _pops, iopt);
942 _expr.template valued_check<value_type>();
943 _prev_omega_K = omega_K;
944 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
945 if (_Xh.get_basis().is_continuous()) {
946 std::set<size_type> ext_dis_idof;
947 _Xh.get_parent_subgeo_owner_dis_indexes (ext_dis_idof);
948 _is_marked.resize (_Xh.ownership());
949 std::fill (_is_marked.begin(), _is_marked.end(),
false);
950 _is_marked.set_dis_indexes (ext_dis_idof);
954template <
class Value>
955typename std::enable_if<
965 Eigen::Matrix<Value,Eigen::Dynamic,1> value;
966 _expr.evaluate (omega_K, K, value);
968 const piola_fem<float_type>& pf = _Xh.get_basis().get_piola_fem();
969 if (pf.transform_need_piola()) {
970 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>&
piola = _pops.get_piola (omega_K, K);
971 for (
size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
973 pf.inv_transform (
piola[loc_inod], value[loc_inod], value[loc_inod]);
976 _Xh.get_basis().compute_dofs (K, value, uk);
979template <
class Value>
980typename std::enable_if<
985 const geo_type& omega_K,
987 vector_element_type& uk)
const
989 using T = scalar_type;
994 switch (_Xh.valued_tag()) {
995#define _RHEOLEF_switch(VALUED,VALUE) \
996 case space_constant::VALUED: \
997 (*this).template evaluate_internal<VALUE> (omega_K, K, uk); break;
1006#undef _RHEOLEF_switch
1007 default:
error_macro (
"unexpected argument valued tag="<<_Xh.valued_tag());
1017 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.
dis_ie()) {
1023 (*this).template evaluate_internal<value_type> (omega_K, K, uk);
1027 if (_Xh.get_basis().is_continuous()) {
1028 std::vector<size_type> dis_idx;
1029 _Xh.dis_idof (K, dis_idx);
1031 "invalid sizes: dis_idx.size="<<dis_idx.size()<<
", uk.size="<<uk.size());
1033 for (
size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; ++loc_idof) {
1034 size_type dis_idof = dis_idx [loc_idof];
1035 size_type iproc = _Xh.get_parent_subgeo_owner (dis_idof);
1036 if (iproc != my_proc || _is_marked.dis_at (dis_idof)) {
1041 _is_marked.dis_entry (dis_idof) =
true;
1045 _prev_omega_K = omega_K;
1046 _prev_K_dis_ie = K.
dis_ie();
1087 { base1::data().evaluate (omega_K, K, uk); }
1097typename std::enable_if<
1108 ,
typename Expr::memory_type>& Xh,
1116template <
class Expr>
1118typename std::enable_if<
1119 details::is_field_function<Expr>::value
1120 ,details::field_lazy_terminal_interpolate<
1121 details::field_expr_v2_nonlinear_terminal_function<Expr>
1138template<
class T,
class M>
1147template <
class T,
class M,
class Expr>
1149typename std::enable_if<
1150 details::is_field_wdof<Expr>::value
1151 && ! details::is_field<Expr>::value
field::size_type size_type
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
const geo_basic< T, M > & level_set() const
rep::memory_type memory_type
rep::scalar_type scalar_type
band_basic< float_type, memory_type > band_type
geo_element::size_type size_type
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
void initialize(const geo_type &omega_K)
typename field_type::space_type space_type
field_lazy_terminal_field(const field_type &uh)
typename field_type::geo_type geo_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
geo_element::size_type size_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename Expr::geo_type geo_type
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
field_lazy_terminal_integrate_band_rep(const band_type &gh, const Expr &expr, const integrate_option &iopt)
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
vector_element_type _prev_uk
typename Expr::space_type space_type
field_lazy_terminal_integrate_band(const band_type &gh, const Expr &expr, const integrate_option &iopt)
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename rep::size_type size_type
typename rep::band_type band_type
typename rep::space_type space_type
typename rep::geo_type geo_type
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
typename rep::vector_element_type vector_element_type
smart_pointer_nocopy< rep > base1
field_lazy_terminal_integrate_band_rep< Expr > rep
typename rep::scalar_type scalar_type
band_basic< float_type, memory_type > band_type
geo_element::size_type size_type
void initialize(const geo_type &omega_K) const
field_lazy_terminal_integrate_rep(const Expr &expr, const integrate_option &iopt)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename Expr::geo_type geo_type
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
disarray< int, memory_type > _is_on_domain
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
const geo_type & get_default_geo() const
typename float_traits< scalar_type >::type float_type
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
vector_element_type _prev_uk
typename Expr::space_type space_type
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename rep::size_type size_type
typename rep::band_type band_type
field_lazy_terminal_integrate_rep< Expr > rep
typename rep::space_type space_type
typename rep::geo_type geo_type
field_lazy_terminal_integrate(const Expr &expr, const integrate_option &iopt)
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
typename rep::vector_element_type vector_element_type
smart_pointer_nocopy< rep > base1
typename rep::scalar_type scalar_type
field_lazy_terminal_integrate(std::string domname, const Expr &expr, const integrate_option &iopt)
field_lazy_terminal_integrate(const geo_type &domain, const Expr &expr, const integrate_option &iopt)
geo_element::size_type size_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
std::enable_if< is_undeterminated< Value >::value, void >::type evaluate_internal(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
disarray< int, memory_type > _is_marked
std::enable_if<!is_undeterminated< Value >::value, void >::type evaluate_internal(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
geo_basic< float_type, memory_type > geo_type
band_type get_band() const
piola_on_pointset< float_type > _pops
const geo_type & get_geo() const
const space_type & get_space() const
band_basic< float_type, memory_type > band_type
typename Expr::value_type value_type
vector_element_type _prev_uk
field_lazy_terminal_interpolate_rep(const space_type &Xh, const Expr &expr)
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
field_lazy_terminal_interpolate_rep< Expr > rep
typename rep::size_type size_type
typename rep::band_type band_type
typename rep::space_type space_type
typename rep::geo_type geo_type
field_lazy_terminal_interpolate(const space_type &Xh, const Expr &expr)
band_type get_band() const
const geo_type & get_geo() const
const space_type & get_space() const
typename rep::vector_element_type vector_element_type
smart_pointer_nocopy< rep > base1
typename rep::scalar_type scalar_type
see the disarray page for the full documentation
geo_basic< float_type, memory_type > geo_type
const geo_type & get_geo() const
const space_type & get_space() const
space_basic< float_type, memory_type > space_type
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
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type n_subgeo(size_type subgeo_dim) const
size_type subgeo_dis_index(size_type subgeo_dim, size_type i) const
see the integrate_option page for the full documentation
see the reference_element page for the full documentation
static const size_type max_side_by_variant
reference_element subgeo(size_type subgeo_dim, size_type loc_isid) const
size_type dimension() const
variant_type variant() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
size_type n_subgeo(size_type subgeo_dim) const
see the tensor3 page for the full documentation
see the tensor4 page for the full documentation
see the tensor page for the full documentation
#define trace_macro(message)
#define error_macro(message)
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_switch(VALUED, VALUE)
void compute_idof_S2idof_K(const basis_basic< T > &b, const reference_element &hat_K, std::array< std::vector< size_t >, reference_element::max_side_by_variant > &idof_S2idof_K, std::array< size_t, reference_element::max_side_by_variant > &ndof_S, size_t &ndof_K)
This file is part of Rheolef.
field_basic< T, M > lazy_interpolate(const space_basic< T, M > &X2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
std::enable_if< details::is_field_expr_quadrature_arg< Expr >::value, details::field_lazy_terminal_integrate< Expr > >::type lazy_integrate(const typename Expr::geo_type &domain, const Expr &expr, const integrate_option &iopt=integrate_option())
see the integrate page for the full documentation
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
helper for std::complex<T>: get basic T type