21#include "rheolef/geo_domain.h"
22#include "rheolef/test.h"
29template <
class T,
class M>
33 _is_inside_on_local_sides(false),
38template <
class T,
class M>
42 _is_inside_on_local_sides (x._is_inside_on_local_sides),
43 _is_on_band (x._is_on_band),
48template <
class T,
class M>
52 trace_macro (
"*** PHYSICAL ASSIGN OF TEST_REP ***");
63template <
class T,
class M>
71 _fops.initialize (get_vf_space().get_basis(), pops);
73template <
class T,
class M>
81 _fops.initialize (get_vf_space().get_basis(), pops);
83template <
class T,
class M>
92 _fops.initialize (get_vf_space().get_basis(), pops);
97template <
class T,
class M,
class Value, details::differentiate_option::type Diff>
98struct evaluate_band_continued {
102 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola,
105 const std::vector<size_t>& dis_inod_K,
106 const std::vector<size_t>& dis_inod_L,
107 const details::differentiate_option& gopt,
108 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
110 fatal_macro (
"on band: unexpected call to differential operator (id="<<Diff<<
")");
113template <
class T,
class M,
class Value>
114struct evaluate_band_continued<
T,
M,Value,details::differentiate_option::none> {
116 const test_rep<T,M>& obj,
117 const geo_basic<T,M>& omega_K,
118 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola,
119 const reference_element& hat_K,
120 const reference_element& tilde_L,
121 const std::vector<size_t>& dis_inod_K,
122 const std::vector<size_t>& dis_inod_L,
123 const details::differentiate_option& gopt,
124 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
126 Eigen::Matrix<Value,Eigen::Dynamic,1> value_i;
127 size_t loc_nnod = value.rows();
128 for (
size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
130 obj._fops.get_basis_on_pointset().get_basis().evaluate (tilde_L, tilde_xi, value_i);
131 value.row(loc_inod) = value_i.transpose();
136template <
class T,
class GradValue>
140 const tensor_basic<T>& invDF,
141 const tensor_basic<T>& P,
142 const details::differentiate_option& gopt,
143 const GradValue& hat_grad_u,
146 fatal_macro(
"band_grad_post: unsupported gradient type: "<<typename_macro(GradValue));
152 const tensor_basic<T>& invDF,
153 const tensor_basic<T>& P,
154 const details::differentiate_option& gopt,
155 const point_basic<T>& tilde_grad_u,
158 grad_u = invDF.trans_mult (tilde_grad_u);
167 const tensor_basic<T>& invDF,
168 const tensor_basic<T>& P,
169 const details::differentiate_option& gopt,
170 const tensor_basic<T>& tilde_grad_u,
173 grad_u = tilde_grad_u*invDF;
174 if (gopt.symmetrized) {
178 if (! gopt.symmetrized) {
225template <
class T,
class M,
class Value>
226struct evaluate_band_continued<
T,
M,Value,details::differentiate_option::gradient> {
228 const test_rep<T,M>& obj,
229 const geo_basic<T,M>& omega_K,
230 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola,
231 const reference_element& hat_K,
232 const reference_element& tilde_L,
233 const std::vector<size_t>& dis_inod_K,
234 const std::vector<size_t>& dis_inod_L,
235 const details::differentiate_option& gopt,
236 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
238 Eigen::Matrix<Value,Eigen::Dynamic,1> value_i;
239 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1> DG;
243 size_t d = omega_K.dimension();
244 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
245 Eigen::Matrix<Value,Eigen::Dynamic,1> tilde_value_i;
246 piola_transformation (obj._gh.level_set(), obj._fops.get_piola_on_pointset().get_basis_on_pointset(), hat_K, dis_inod_K, x);
247 size_t loc_nnod = value.rows();
248 size_t loc_ndof = value.cols();
249 for (
size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
253 obj._fops.get_basis_on_pointset().get_basis().grad_evaluate (tilde_L, tilde_xi, tilde_value_i);
257 jacobian_piola_transformation (obj._gh.band(), obj._fops.get_piola_on_pointset().get_basis_on_pointset().get_basis(), tilde_L, dis_inod_L, tilde_xi, DF);
264 for (
size_t loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
265 const Value& hat_grad_u = tilde_value_i [loc_jdof];
266 Value&
grad_u = value (loc_inod, loc_jdof);
267 band_grad_post (invDF, P, gopt, hat_grad_u,
grad_u);
278template <
class T,
class M>
281 const geo_basic<T,M>& space_geo,
282 const geo_basic<T,M>& omega_K,
283 const geo_element& K_in);
285template <
class T,
class M>
286template <
class Value, details::differentiate_option::type Diff>
292 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
306 if (!use_bgd2dom && !use_dom2bgd) {
309 1 + K.
dimension() == omega_K.get_background_geo().map_dimension() &&
310 _V.get_basis().is_discontinuous();
311 if (!use_dg_on_sides) {
313 _fops.template evaluate<M,Value,Diff> (omega_K, K, gopt, value);
320 check_macro (L_map_d == omega_L_eff.map_dimension(),
321 "unexpected dimension for side S="<<K.
name()<<K.
dis_ie()<<
" in domain "<<omega_L_eff.name());
323 "unexpected non-boundary side S="<<K.
name()<<K.
dis_ie()<<
" in domain "<<omega_L_eff.name());
325 check_macro (L_dis_ie != std::numeric_limits<size_type>::max(),
326 "unexpected isolated side S="<<K.
name()<<K.
dis_ie());
327 const geo_element& L_eff = omega_L_eff.dis_get_geo_element (L_map_d, L_dis_ie);
330 _fops.template evaluate_on_side<M,Value,Diff> (omega_L_eff, L_eff, sid, gopt, value);
332 }
else if (use_dom2bgd) {
335 1 + K.
dimension() == _V.get_geo().get_background_geo().map_dimension() &&
336 _V.get_basis().is_discontinuous();
337 if (!use_dg_on_sides) {
339 _fops.template evaluate<M,Value,Diff> (omega_K, K, gopt, value);
345 check_macro (L_space_map_d == _V.get_geo().map_dimension(),
346 "unexpected dimension for side S="<<K_space.
name()<<K_space.
dis_ie()<<
" in domain "<<_V.get_geo().name());
348 "unexpected non-boundary side S="<<K_space.
name()<<K_space.
dis_ie()<<
" in domain "<<_V.get_geo().name());
350 check_macro (L_space_dis_ie != std::numeric_limits<size_type>::max(),
351 "unexpected isolated side S="<<K_space.
name()<<K_space.
dis_ie());
352 const geo_element& L_space = _V.get_geo().dis_get_geo_element (L_space_map_d, L_space_dis_ie);
355 _fops.template evaluate_on_side<M,Value,Diff> (_V.get_geo(), L_space, sid, gopt, value);
368 omega_K.get_background_geo().name() == _V.get_geo().name() ||
369 omega_K.get_background_geo().name() == _V.get_geo().get_background_geo().name(),
370 "unexpected integration domain \"" << omega_K.name()
371 <<
"\" based on geo \"" << omega_K.get_background_geo().name()
372 <<
"\" together with a space on geo \"" << _V.get_geo().name()
373 <<
"\" based on geo \"" << _V.get_geo().get_background_geo().name()<<
"\"");
375 const geo_element& K_eff = omega_K_eff.bgd2dom_geo_element (K);
376 _fops.template evaluate<M,Value,Diff> (omega_K_eff, K_eff, gopt, value);
381 std::vector<size_type> dis_inod_L;
382 std::vector<size_type> dis_inod_K;
383 size_type first_dis_ie = _gh.level_set().sizes().ownership_by_dimension [K.
dimension()].first_index();
386 size_type L_ie = _gh.sid_ie2bnd_ie (K_ie);
392 const piola_fem<T>& pf = _fops.get_basis_on_pointset().get_basis().get_piola_fem();
393 check_macro (! pf.transform_need_piola(),
"band: unsupported piola-based basis");
394 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&
piola = _fops.get_piola_on_pointset().get_piola (_gh.level_set(), K);
396 size_type loc_ndof = get_vf_space().get_basis().ndof (tilde_L);
397 value.resize (loc_nnod, loc_ndof);
398 _gh.level_set().dis_inod (K, dis_inod_K);
399 _gh.band() .dis_inod (L, dis_inod_L);
400 evaluate_band_continued<T,M,Value,Diff> eval_band_cont;
401 eval_band_cont (*
this, omega_K,
piola, hat_K, tilde_L, dis_inod_K, dis_inod_L, gopt, value);
407template <
class T,
class M>
408template <
class Value, details::differentiate_option::type Diff>
415 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
417 size_type space_map_d = get_vf_space().get_geo().map_dimension();
419 bool is_broken = get_vf_space().get_geo().is_broken() && space_map_d + 1 == K.
dimension();
422 _fops.template evaluate_on_side<M,Value,Diff> (omega_K, K, sid, gopt, value);
424 check_macro(_is_inside_on_local_sides,
"unexpected broken case");
427 const geo_element& S = omega_K.dis_get_geo_element (space_map_d, dis_isid);
429 const geo_element* dom_S_ptr = (!have_bgd_dom) ? &S : &(get_vf_space().get_geo().bgd2dom_geo_element(S));
431 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> value_sid;
432 _fops.template evaluate<M,Value,Diff> (omega_K, dom_S, gopt, value_sid);
441 size_type loc_jsid_ndof = get_vf_space().get_basis().ndof (hat_Sj);
443 start_loc_isid_idof = loc_ndof;
444 loc_isid_ndof = loc_jsid_ndof;
446 loc_ndof += loc_jsid_ndof;
448 value.resize (loc_nnod, loc_ndof);
449 value.fill (Value());
450 for (
size_type loc_isid_jdof = 0; loc_isid_jdof < loc_isid_ndof; ++loc_isid_jdof) {
451 size_type loc_jdof = start_loc_isid_idof + loc_isid_jdof;
452 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
453 value(loc_inod,loc_jdof) = value_sid (loc_inod, loc_isid_jdof);
458 check_macro (!_is_on_band,
"evaluate_on_side: not yet on band");
464template <
class T,
class M>
465template <
class Value>
472 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
473 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
474 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
476 check_macro (value0.rows() == value1.rows(),
"invalid node number");
480 value.resize (loc_nnod, loc_ndof0 + loc_ndof1);
481 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
482 for (
size_type loc_idof = 0; loc_idof < loc_ndof0; ++loc_idof) {
483 value (loc_inod,loc_idof) = value0 (loc_inod,loc_idof);
485 for (
size_type loc_idof = 0; loc_idof < loc_ndof1; ++loc_idof) {
486 value (loc_inod,loc_ndof0+loc_idof) = value1 (loc_inod,loc_idof);
494#define _RHEOLEF_instanciation_value_diff(T,M,Value,Diff) \
495template void test_rep<T,M>::evaluate<Value,Diff> ( \
496 const geo_basic<T,M>& omega_K, \
497 const geo_element& K, \
498 const details::differentiate_option& gopt, \
499 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const; \
500template void test_rep<T,M>::evaluate_on_side<Value,Diff> ( \
501 const geo_basic<T,M>& omega_K, \
502 const geo_element& K, \
503 const side_information_type& sid, \
504 const details::differentiate_option& gopt, \
505 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const; \
507#define _RHEOLEF_instanciation_value(T,M,Value) \
508template void test_rep<T,M>::local_dg_merge_on_side ( \
509 const geo_basic<T,M>& omega_K, \
510 const geo_element& S, \
511 const geo_element& K0, \
512 const geo_element& K1, \
513 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0, \
514 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1, \
515 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const; \
516_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::none) \
517_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::gradient) \
518_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::divergence) \
519_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::curl) \
521#define _RHEOLEF_instanciation(T,M) \
522template class test_rep<T,M>; \
523_RHEOLEF_instanciation_value(T,M,T) \
524_RHEOLEF_instanciation_value(T,M,point_basic<T>) \
525_RHEOLEF_instanciation_value(T,M,tensor_basic<T>) \
526_RHEOLEF_instanciation_value(T,M,tensor3_basic<T>) \
527_RHEOLEF_instanciation_value(T,M,tensor4_basic<T>) \
530#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciation(T, M, A)
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the Float page for the full documentation
abstract base interface class
virtual size_type variant() const =0
virtual std::string name() const =0
generic mesh with rerefence counting
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
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
reference_element side(size_type loc_isid) const
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
band_basic< float_type, M > _gh
fem_on_pointset< float_type > _fops
bool _is_inside_on_local_sides
test_rep(const space_type &V)
void evaluate_on_side(const geo_basic< T, M > &omega_K, const geo_element &K, const side_information_type &sid, const details::differentiate_option &gopt, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
void local_dg_merge_on_side(const geo_basic< T, M > &omega_K, const geo_element &S, const geo_element &K0, const geo_element &K1, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value0, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value1, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
void evaluate(const geo_basic< T, M > &omega_K, const geo_element &K, const details::differentiate_option &gopt, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
test_rep< T, M > & operator=(const test_rep< T, M > &)
#define trace_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)")
This file is part of Rheolef.
void piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
tensor_basic< T > pseudo_inverse_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
void map_projector(const tensor_basic< T > &DF, size_t d, size_t map_d, tensor_basic< T > &P)
void jacobian_piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &DF)
point_basic< T > inverse_piola_transformation(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod, const point_basic< T > &x)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
const geo_element & assembly2space_geo_element(const geo_basic< T, M > &space_geo, const geo_basic< T, M > &omega_K, const geo_element &K_in)