field_vf_assembly.h
Go to the documentation of this file.
#ifndef _RHEO_FIELD_VF_ASSEMBLY_H
#define _RHEO_FIELD_VF_ASSEMBLY_H
#include "rheolef/field.h"
#include "rheolef/test.h"
#include "rheolef/integrate_option.h"
#include "rheolef/field_expr_quadrature.h"
/*
let:
l(v) = int_omega expr(v) dx
The integrals are evaluated over each element K of the domain
by using a quadrature formulae given by iopt
expr(v) is a linear expression with respect to the test-function v
The test-function v is replaced by each of the basis function of
the corresponding finite element space Xh: (phi_i), i=0..dim(Xh)-1
The integrals over K are transformed on the reference element with
the piola transformation:
F : hat_K ---> K
hat_x |--> x = F(hat_x)
exemples:
1) expr(v) = v
int_K phi_i(x) dx
= int_{hat_K} hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
= sum_q hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
The value(q,i) = (hat_phi_i(hat_xq)) are evaluated on time for all over the
reference element hat_K and for the given quadrature formulae by:
expr.initialize (omega, quad);
This expression is represented by the 'test' class (see test.h)
2) expr(v) = f*v
int_K f(x)*phi_i(x) dx
= int_{hat_K} f(F(hat_x)*hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
= sum_q f(F(hat_xq))*hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
On each K, the computation needs two imbricated loops in (q,i).
The expresion is represented by a tree at compile-time.
The first subexpr 'f' is represented by field_vf_expr_terminal<f> : it evaluates :
value1(q,i) = f(F(hat_xq)) : the values depends only of q and are independent of i.
They could be evaluated juste before the 'i' loop.
As the field_vf_expr_terminal<> is general and can handle also more complex expressions
involving fields with various approximations, all the values on K are evaluated
in a specific 'q' loop by
subexpr1.element_initialize (K);
The second subexpr 'phi_i' is represdented by the 'test' class (see before).
value2(q,i) = (hat_phi_i(hat_xq))
The '*' performs on the fly the product
value(q,i) = value1(q,i)*value2(q,i)
3) expr(v) = dot(f,grad(v)) dx
The 'f' function is here vector-valued.
The expresion is represented by a tree at compile-time :
dot -> f
-> grad -> v
The 'grad' node returns
value(q,i) = trans(inv(DF(hat_wq))*grad_phi_i(hat_xq) that is vector-valued
The grad_phi values are obtained by a grad_value(q,i) method on the 'test' class.
The 'dot' performs on the fly the product
value(q,i) = ot (value1(q,i), value2(q,i))
This approch generilize for an expression tree.
*/
namespace rheolef {
// ====================================================================
// common implementation for integration on a band or an usual domain
// ====================================================================
template <class T, class M>
template <class Expr>
void
field_basic<T,M>::assembly_internal (
const geo_basic<T,M>& omega,
const geo_basic<T,M>& band,
const band_basic<T,M>& gh,
const Expr& expr,
const integrate_option& iopt,
bool is_on_band)
{
// ------------------------------------
// 0) initialize the loop
// ------------------------------------
const space_basic<T,M>& Xh = expr.get_vf_space();
const geo_basic<T,M>& bgd_omega = Xh.get_constitution().get_background_geo();
if (is_on_band) {
check_macro (band.get_background_geo() == bgd_omega,
"assembly: incompatible integration domain "<<omega.name() << " and test function based domain "
<< bgd_omega.name());
}
resize (Xh, 0);
if (is_on_band) {
expr.initialize (gh, iopt);
} else {
expr.initialize (omega, iopt);
}
expr.template valued_check<T>();
vec<T,M>& u = set_u();
vec<T,M>& b = set_b();
std::vector<size_type> dis_idx;
Eigen::Matrix<T,Eigen::Dynamic,1> lk;
size_type d = omega.dimension();
size_type map_d = omega.map_dimension();
if (Xh.get_constitution().is_discontinuous()) Xh.get_constitution().neighbour_guard();
// ------------------------------------
// 1) loop on elements
// ------------------------------------
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 dofs
// ------------------------------------
if (!is_on_band) {
Xh.get_constitution().assembly_dis_idof (omega, K, dis_idx);
} else { // on a band
size_type L_ie = gh.sid_ie2bnd_ie (ie);
const geo_element& L = band [L_ie];
Xh.dis_idof (L, dis_idx);
}
// ------------------------------------
// 1.2) locally compute values
// ------------------------------------
expr.evaluate (omega, K, lk);
// -----------------------------------------
// 1.3) assembly local lk in global field lh
// -----------------------------------------
check_macro (dis_idx.size() == size_type(lk.size()),
"incompatible sizes dis_idx("<<dis_idx.size()<<") and lk("<<lk.size()<<")");
for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
const T& value = lk [loc_idof];
size_type dis_idof = dis_idx [loc_idof];
size_type dis_iub = Xh.dis_idof2dis_iub(dis_idof);
if (Xh.dis_is_blocked(dis_idof)) b.dis_entry(dis_iub) += value;
else u.dis_entry(dis_iub) += value;
}
}
// -----------------------------------------
// 2) finalize the assembly process
// -----------------------------------------
u.dis_entry_assembly(set_add_op<T,T>());
b.dis_entry_assembly(set_add_op<T,T>());
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
const geo_basic<T,M>& omega,
const Expr& expr,
const integrate_option& iopt)
{
assembly_internal (omega, omega, band_basic<T,M>(), expr, iopt, false);
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
const band_basic<T,M>& gh,
const Expr& expr,
const integrate_option& iopt)
{
assembly_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
}
}// namespace rheolef
#endif // _RHEO_FIELD_VF_ASSEMBLY_H
rheolef::geo_basic< T, M >
rheolef::field_basic::assembly
void assembly(const geo_basic< T, M > &dom, const Expr &expr, const integrate_option &iopt)
Definition: field_vf_assembly.h:170
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
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::value
rheolef::std value
rheolef::space_basic
the finite element space
Definition: space.h:352
band
see the band page for the full documentation
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::band_basic
Definition: band.h:95
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
u
Definition: leveque.h:25
mkgeo_obstacle.L
L
Definition: mkgeo_obstacle.sh:133
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
rheolef::field_basic::assembly_internal
void assembly_internal(const geo_basic< T, M > &dom, const geo_basic< T, M > &band, const band_basic< T, M > &gh, const Expr &expr, const integrate_option &qopt, bool is_on_band)
Definition: field_vf_assembly.h:95
rheolef::set_add_op
Definition: msg_util.h:61
rheolef::field_basic< T, rheo_default_memory_model >::size_type
std::size_t size_type
Definition: field.h:239
T
Expr1::float_type T
Definition: field_expr.h:218