characteristic.cc
Go to the documentation of this file.
#include "rheolef/characteristic.h"
#include "rheolef/fem_on_pointset.h"
#include "rheolef/rheostream.h"
#include "rheolef/band.h"
#include "rheolef/field_evaluate.h"
// TODO: replace piola_on_quad aka basis_on_pointset by piola_on_poinset
// ===========================================================================
// allocator
// ===========================================================================
namespace rheolef {
template<class T, class M>
characteristic_rep<T,M>::characteristic_rep(const field_basic<T,M>& dh)
: _dh(dh),
_coq_map()
{
// automatically add the "boundary" domain, if not already present:
_dh.get_geo().boundary();
}
} // namespace rheolef
// ===========================================================================
// externs
// ===========================================================================
namespace rheolef { namespace details {
template<class T, class M>
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, // x has been already localized in K
disarray<index_set,M>& ie2dis_ix, // K -> list of ix
disarray<point_basic<T>,M>& hat_y);
// ===========================================================================
// integrate(omega,compose (uh,X)) returns a field:
// It is done in two passes:
// - the first one is symbolic: localize the moved quadrature points yq=xq+X.dh(xq)
// this is done one time for all
// - the second pass is valued: compte uh(yq) during integration
// ===========================================================================
template <class T, class M>
void
integrate_pass1_symbolic (
// input :
const space_basic<T,M>& Xh,
const field_basic<T,M>& dh,
const piola_on_pointset<T>& pops,
// output :
disarray<index_set,M>& ie2dis_ix,
disarray<point_basic<T>,M>& hat_y,
disarray<point_basic<T>,M >& yq)
{
// ----------------------------------------------------------------------
// 1) set the quadrature formulae
// ----------------------------------------------------------------------
typedef typename space_basic<T,M>::size_type size_type;
const geo_basic<T,M>& omega = Xh.get_geo();
check_macro (omega == dh.get_geo(), "characteristic: unsupported different meshes for flow ("
<< dh.get_geo().name() << ") and convected field (" << omega.name() << ")");
if (omega.order() > 1) {
warning_macro ("characteristic: high-order curved elements not yet supported (treated as first order)");
}
fem_on_pointset<T> d_fops;
d_fops.initialize (dh.get_space().get_constitution().get_basis(), pops);
// ----------------------------------------------------------------------
// 2) size of the the disarray of all quadrature nodes
// ----------------------------------------------------------------------
// NOTE: since these nodes are moved (convected) by the flow associated to the characteristic
// they can change from partition and go to an element managed by another proc
// thus, in order to group comms, a whole disarray of quadrature point is allocated:
// In sequential mode (or with only one proc), this can do in an element by element way
// with less memory.
size_type dim = omega.dimension();
size_type map_dim = omega.map_dimension();
size_type sz = 0;
size_type dis_sz = 0;
const basis_on_pointset<T>& piola_on_quad = pops.get_basis_on_pointset();
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
distributor ige_ownership = omega.sizes().ownership_by_variant [variant];
if (ige_ownership.dis_size() == 0) continue;
reference_element hat_K (variant);
size_type nq = piola_on_quad.nnod (hat_K);
sz += nq*ige_ownership.size();
dis_sz += nq*ige_ownership.dis_size();
}
distributor xq_ownership (dis_sz, omega.comm(), sz);
// ----------------------------------------------------------------------
// 3) build the disarray of all quadrature nodes; xq
// ----------------------------------------------------------------------
disarray<point_basic<T>,M> xq (xq_ownership);
std::vector<size_type> dis_inod;
for (size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega[ie];
const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
for (size_type q = 0, nq = piola.size(); q < nq; q++, ixq++) {
xq[ixq] = piola[q].F;
}
}
// ----------------------------------------------------------------------
// 4) interpolate dh on the xq nodes: dq(i)=dh(xq(i))
// ----------------------------------------------------------------------
dh.dis_dof_update();
disarray<point_basic<T>,M> dq (xq_ownership);
disarray<size_type,M> ix2dis_ie (xq_ownership);
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> d_value;
for (size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega[ie];
field_evaluate (dh, d_fops, omega, K, d_value);
for (size_type q = 0, nq = d_value.size(); q < nq; q++, ixq++) {
dq[ixq] = d_value [q];
ix2dis_ie[ixq] = ie; // will be a guest for locate(yq)
}
}
// ----------------------------------------------------------------------
// 5) build the disarray of all moved quadrature nodes yq(i) = xq(i)+dq(i)
// with the constraint to remain inside Omega
// ----------------------------------------------------------------------
yq.resize (xq_ownership);
omega.trace_move (xq, dq, ix2dis_ie, yq);
// ----------------------------------------------------------------------
// 6.a) symbolic interpolation-like pass: localize the yq nodes
// ----------------------------------------------------------------------
interpolate_pass1_symbolic (omega, yq, ix2dis_ie, ie2dis_ix, hat_y);
}
} // namespace details
// ========================================================================
// symbolic pass 1 is stored one time for all in the characteristic class
// ========================================================================
template<class T, class M>
const characteristic_on_quadrature<T,M>&
characteristic_rep<T,M>::get_pre_computed (
const space_basic<T,M>& Xh,
const field_basic<T,M>& dh,
const piola_on_pointset<T>& pops) const
{
// get the quadrature used for integration
check_macro (pops.has_quadrature(), "unexpected point set, expect quadrature point set");
const quadrature<T>& quad = pops.get_quadrature();
const quadrature_option& qopt = quad.get_options();
check_macro (qopt.get_family() == quadrature_option::gauss_lobatto,
"integrate characteristics HINT: use Gauss-Lobatto quadrature)");
check_macro (qopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
"integrate characteristics HINT: invalid unset order");
// search if already used & memorized
std::string label = qopt.get_family_name() + "(" + itos(qopt.get_order()) + ")";
typename map_type::const_iterator iter = _coq_map.find (label);
if (iter != _coq_map.end()) {
return (*iter).second;
}
// build a new one & memorize it
// search if already used & memorized
characteristic_on_quadrature<T,M> coq (qopt);
characteristic_on_quadrature_rep<T,M>& coq_r = coq.data();
coq_r._qopt = qopt;
coq_r._quad = quad;
coq_r._piola_on_quad = pops.get_basis_on_pointset();
details::integrate_pass1_symbolic (Xh, dh, pops,
coq_r._ie2dis_ix,
coq_r._hat_y,
coq_r._yq);
// output :
std::pair <typename map_type::iterator,bool> iter_b = _coq_map.insert (std::make_pair(label,coq));
typename map_type::iterator iter2 = iter_b.first;
return (*iter2).second;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template class characteristic_rep<T,M>; \
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
}// namespace rheolef
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::geo_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
rheolef::point_basic
Definition: point.h:87
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::details::interpolate_pass1_symbolic
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)
Definition: interpolate.cc:44
rheolef::smart_pointer_base< characteristic_on_quadrature_rep< T, M >, details::constructor_copy< characteristic_on_quadrature_rep< T, M > > >::data
const characteristic_on_quadrature_rep< T, M > & data() const
Definition: smart_pointer.h:266
rheolef::piola::F
point_basic< T > F
Definition: piola.h:79
rheolef::characteristic_on_quadrature_rep::_hat_y
disarray< point_basic< T >, M > _hat_y
Definition: characteristic.h:146
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::characteristic_on_quadrature_rep
Definition: characteristic.h:118
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::piola_on_pointset::get_quadrature
const quadrature< T > & get_quadrature() const
Definition: piola_on_pointset.h:224
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::details::integrate_pass1_symbolic
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)
Definition: characteristic.cc:67
rheolef::quadrature::get_options
const quadrature_option & get_options() const
Definition: quadrature_rep.cc:264
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::characteristic_on_quadrature_rep::_piola_on_quad
basis_on_pointset< T > _piola_on_quad
Definition: characteristic.h:144
rheolef::characteristic_rep::characteristic_rep
characteristic_rep(const field_basic< T, M > &dh)
Definition: characteristic.cc:35
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::characteristic_on_quadrature_rep::_yq
disarray< point_basic< T >, M > _yq
Definition: characteristic.h:147
rheolef::fem_on_pointset::initialize
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
Definition: fem_on_pointset.h:130
rheolef::field_basic
Definition: field_expr_utilities.h:38
rheolef::characteristic_rep::get_pre_computed
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
Definition: characteristic.cc:157
rheolef::characteristic_on_quadrature_rep::_qopt
quadrature_option _qopt
Definition: characteristic.h:142
rheolef::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::distributor::dis_size
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
rheolef::fem_on_pointset
Definition: fem_on_pointset.h:117
rheolef::characteristic_on_quadrature
Definition: characteristic.h:150
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
rheolef::characteristic_rep::_dh
field_basic< T, M > _dh
Definition: characteristic.h:187
rheolef::characteristic_on_quadrature_rep::_quad
quadrature< T > _quad
Definition: characteristic.h:143
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::basis_on_pointset
Definition: basis_on_pointset.h:191
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::basis_on_pointset::nnod
size_type nnod(reference_element hat_K) const
Definition: basis_on_pointset.h:263
rheolef::field_evaluate
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)
Definition: field_evaluate.cc:39
rheolef::characteristic_on_quadrature_rep::_ie2dis_ix
disarray< index_set, M > _ie2dis_ix
Definition: characteristic.h:145
size_type
field::size_type size_type
Definition: branch.cc:425
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::distributed
distributed
Definition: asr.cc:228
mkgeo_ball.dim
dim
Definition: mkgeo_ball.sh:307
rheolef::piola_on_pointset::has_quadrature
bool has_quadrature() const
Definition: piola_on_pointset.h:217
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::quadrature
Definition: quadrature.h:185
M
Expr1::memory_type M
Definition: vec_expr_v2.h:385
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::piola_on_pointset::get_basis_on_pointset
const basis_on_pointset< T > & get_basis_on_pointset() const
Definition: piola_on_pointset.h:210
rheolef::piola_on_pointset::get_piola
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
Definition: piola_on_pointset.h:239
rheolef::piola
Definition: piola.h:67