1 #ifndef _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
2 #define _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
3 #include "rheolef/field.h"
24 #include "rheolef/test.h"
25 #include "rheolef/quadrature.h"
51 namespace rheolef {
namespace details {
53 template <
class T,
class M,
class Expr,
class Result>
67 "integrate: the quadrature formulae order may be chosen (HINT: use qopt.set_order(n))");
71 pops.
initialize (omega.get_piola_basis(), quad, iopt);
72 expr.initialize (pops, iopt);
73 expr.template valued_check<Result>();
75 Eigen::Matrix<Result,Eigen::Dynamic,1>
value;
80 for (
size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
81 const geo_element& K = omega.get_geo_element (map_d, ie);
88 expr.evaluate (omega, K,
value);
89 const Eigen::Matrix<T,Eigen::Dynamic,1>& w = pops.
get_weight (omega, K);
91 "incompatible sizes w("<<w.size()<<
") and value("<<
value.size()<<
")");
95 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
96 result +=
value[loc_inod]*w[loc_inod];
99 #ifdef _RHEOLEF_HAVE_MPI
101 result = mpi::all_reduce (omega.geo_element_ownership(0).comm(), result, std::plus<Result>());
103 #endif // _RHEOLEF_HAVE_MPI
107 #endif // _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H