Loading...
Searching...
No Matches
field_expr_value_assembly.h
Go to the documentation of this file.
#ifndef _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
#define _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
#include "rheolef/field.h"
#include "rheolef/test.h"
#include "rheolef/quadrature.h"
/*
let:
I = int_omega expr(x) dx
The integrals are evaluated over each element K of the domain omega
by using a quadrature formulae given by iopt
expr(x) is a nonlinear expression
The integrals over K are transformed on the reference element with
the piola transformation:
F : hat_K ---> K
hat_x |--> x = F(hat_x)
int_K expr(x) dx
= int_{hat_K} expr(F(hat_x) det(DF(hat_x)) d hat_x
= sum_q expr(F(hat_xq)) det(DF(hat_xq)) hat_wq
On each K, the computation needs a loop in q.
The expresion is represented by a tree at compile-time.
The 'expr' is represented by field_expr_terminal<f> : it evaluation gives :
value1(q) = expr(F(hat_xq))
This approch generilizes to an expression tree.
*/
namespace rheolef { namespace details {
template <class T, class M, class Expr, class Result>
void
field_expr_v2_value_assembly (
const geo_basic<T,M>& omega,
const Expr& expr0,
const integrate_option& iopt,
Result& result)
{
Expr expr = expr0; // so could expr.initialze()
typedef typename geo_basic<T,M>::size_type size_type;
// ------------------------------------
// 0) initialize the loop
// ------------------------------------
quadrature<T> quad;
check_macro (iopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
"integrate: the quadrature formulae order may be chosen (HINT: use qopt.set_order(n))");
quad.set_order (iopt.get_order());
quad.set_family (iopt.get_family());
piola_on_pointset<T> pops;
pops.initialize (omega.get_piola_basis(), quad, iopt);
expr.initialize (pops, iopt);
expr.template valued_check<Result>();
size_type map_d = omega.map_dimension();
Eigen::Matrix<Result,Eigen::Dynamic,1> value;
// ------------------------------------
// 1) loop on elements
// ------------------------------------
result = Result(0);
for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
const geo_element& K = omega.get_geo_element (map_d, ie);
// ------------------------------------
// 1.1) locally compute values
// ------------------------------------
// locally evaluate int_K expr dx
// with loop on a quadrature formulae
// and accumulate in result
expr.evaluate (omega, K, value);
const Eigen::Matrix<T,Eigen::Dynamic,1>& w = pops.get_weight (omega, K);
check_macro (w.size() == value.size(),
"incompatible sizes w("<<w.size()<<") and value("<<value.size()<<")");
// TODO: DVT_BLAS1 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic> ri = value.transpose()*w;
// => compil pb with Eigen when Result is point, tensor, etc:
size_type loc_nnod = value.size();
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
result += value[loc_inod]*w[loc_inod];
}
}
#ifdef _RHEOLEF_HAVE_MPI
if (is_distributed<M>::value) {
result = mpi::all_reduce (omega.geo_element_ownership(0).comm(), result, std::plus<Result>());
}
#endif // _RHEOLEF_HAVE_MPI
}
}}// namespace rheolef::details
#endif // _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition integrate_option.h:125
rheolef::integrate_option::get_family
family_type get_family() const
Definition integrate_option.h:248
rheolef::integrate_option::get_order
size_t get_order() const
Definition integrate_option.h:242
rheolef::piola_on_pointset
Definition piola_on_pointset.h:142
rheolef::piola_on_pointset::initialize
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
Definition piola_on_pointset.h:190
rheolef::piola_on_pointset::get_weight
const Eigen::Matrix< T, Eigen::Dynamic, 1 > & get_weight(const geo_basic< T, M > &omega, const geo_element &K) const
Definition piola_on_pointset.h:247
rheolef::quadrature
Definition quadrature.h:186
rheolef::quadrature::set_family
void set_family(family_type ft)
Definition quadrature_rep.cc:292
rheolef::quadrature::set_order
void set_order(size_type order)
Definition quadrature_rep.cc:284
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
Definition compiler_eigen.h:39
rheolef::details::field_expr_v2_value_assembly
void field_expr_v2_value_assembly(const geo_basic< T, M > &omega, const Expr &expr0, const integrate_option &iopt, Result &result)
Definition field_expr_value_assembly.h:55
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::valued_check
bool valued_check() const
Definition field_expr_recursive.h:956
rheolef::is_distributed::value
static const bool value
Definition communicator.h:71