21#include "rheolef/characteristic.h"
22#include "rheolef/fem_on_pointset.h"
23#include "rheolef/rheostream.h"
24#include "rheolef/band.h"
25#include "rheolef/field_evaluate.h"
34template<
class T,
class M>
47namespace rheolef {
namespace details {
49template<
class T,
class M>
52 const geo_basic<T,M>& omega,
53 const disarray<point_basic<T>,
M>& x,
54 const disarray<geo_element::size_type,M>& ix2dis_ie,
55 disarray<index_set,M>& ie2dis_ix,
56 disarray<point_basic<T>,
M>& hat_y);
65template <
class T,
class M>
82 check_macro (omega == dh.
get_geo(),
"characteristic: unsupported different meshes for flow ("
83 << dh.
get_geo().name() <<
") and convected field (" << omega.name() <<
")");
84 if (omega.order() > 1) {
85 warning_macro (
"characteristic: high-order curved elements not yet supported (treated as first order)");
98 size_type map_dim = omega.map_dimension();
104 distributor ige_ownership = omega.sizes().ownership_by_variant [variant];
105 if (ige_ownership.
dis_size() == 0)
continue;
108 sz += nq*ige_ownership.
size();
109 dis_sz += nq*ige_ownership.
dis_size();
111 distributor xq_ownership (dis_sz, omega.comm(), sz);
116 std::vector<size_type> dis_inod;
117 for (
size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
119 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&
piola = pops.
get_piola (omega, K);
130 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> d_value;
131 for (
size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
134 for (
size_type q = 0, nq = d_value.size(); q < nq; q++, ixq++) {
135 dq[ixq] = d_value [q];
143 yq.resize (xq_ownership);
144 omega.trace_move (xq, dq, ix2dis_ie, yq);
155template<
class T,
class M>
156const characteristic_on_quadrature<T,M>&
165 const quadrature_option& qopt = quad.
get_options();
167 check_macro (qopt.get_family() == quadrature_option::gauss_lobatto,
168 "integrate characteristics HINT: use Gauss-Lobatto quadrature)");
169 check_macro (qopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
170 "integrate characteristics HINT: invalid unset order");
173 std::string label = qopt.get_family_name() +
"(" + std::to_string(qopt.get_order()) +
")";
174 typename map_type::const_iterator iter = _coq_map.find (label);
175 if (iter != _coq_map.end()) {
176 return (*iter).second;
191 std::pair <typename map_type::iterator,bool> iter_b = _coq_map.insert (std::make_pair(label,coq));
192 typename map_type::iterator iter2 = iter_b.first;
193 return (*iter2).second;
198#define _RHEOLEF_instanciation(T,M) \
199template class characteristic_rep<T,M>; \
202#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciation(T, M, A)
field::size_type size_type
see the Float page for the full documentation
size_type nnod(reference_element hat_K) const
characteristic_rep(const field_basic< T, M > &dh)
const characteristic_on_quadrature< T, M > & get_pre_computed(const space_basic< T, M > &Xh, const field_basic< T, M > &dh, const piola_on_pointset< T > &pops) const
see the disarray page for the full documentation
see the distributor page for the full documentation
size_type dis_size() const
global and local sizes
size_type size(size_type iproc) const
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
void dis_dof_update(const SetOp &=SetOp()) const
const geo_type & get_geo() const
const space_type & get_space() const
generic mesh with rerefence counting
see the geo_element page for the full documentation
bool has_quadrature() const
const quadrature< T > & get_quadrature() const
const basis_on_pointset< T > & get_basis_on_pointset() const
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
const quadrature_option & get_options() const
see the reference_element page for the full documentation
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
#define warning_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void interpolate_pass1_symbolic(const geo_basic< T, M > &omega, const disarray< point_basic< T >, M > &x, const disarray< geo_element::size_type, M > &ix2dis_ie, disarray< index_set, M > &ie2dis_ix, disarray< point_basic< T >, M > &hat_y)
void integrate_pass1_symbolic(const space_basic< T, M > &Xh, const field_basic< T, M > &dh, const piola_on_pointset< T > &pops, disarray< index_set, M > &ie2dis_ix, disarray< point_basic< T >, M > &hat_y, disarray< point_basic< T >, M > &yq)
This file is part of Rheolef.
void field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
disarray< index_set, M > _ie2dis_ix
disarray< point_basic< T >, M > _hat_y
basis_on_pointset< T > _piola_on_quad
disarray< point_basic< T >, M > _yq