Rheolef  7.1
an efficient C++ finite element environment
field_expr_value_assembly.h
Go to the documentation of this file.
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"
26 /*
27 let:
28 I = int_omega expr(x) dx
29 
30 The integrals are evaluated over each element K of the domain omega
31 by using a quadrature formulae given by iopt
32 
33 expr(x) is a nonlinear expression
34 
35 The integrals over K are transformed on the reference element with
36 the piola transformation:
37 F : hat_K ---> K
38  hat_x |--> x = F(hat_x)
39 
40 int_K expr(x) dx
41 = int_{hat_K} expr(F(hat_x) det(DF(hat_x)) d hat_x
42 = sum_q expr(F(hat_xq)) det(DF(hat_xq)) hat_wq
43 
44 On each K, the computation needs a loop in q.
45 The expresion is represented by a tree at compile-time.
46 The 'expr' is represented by field_expr_terminal<f> : it evaluation gives :
47 value1(q) = expr(F(hat_xq))
48 
49 This approch generilizes to an expression tree.
50 */
51 namespace rheolef { namespace details {
52 
53 template <class T, class M, class Expr, class Result>
54 void
56 const geo_basic<T,M>& omega,
57 const Expr& expr,
58 const integrate_option& iopt,
59 Result& result)
60 {
61  typedef typename geo_basic<T,M>::size_type size_type;
62  // ------------------------------------
63  // 0) initialize the loop
64  // ------------------------------------
65  quadrature<T> quad;
66  check_macro (iopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
67  "integrate: the quadrature formulae order may be chosen (HINT: use qopt.set_order(n))");
68  quad.set_order (iopt.get_order());
69  quad.set_family (iopt.get_family());
71  pops.initialize (omega.get_piola_basis(), quad, iopt);
72  expr.initialize (pops, iopt);
73  expr.template valued_check<Result>();
74  size_type map_d = omega.map_dimension();
75  Eigen::Matrix<Result,Eigen::Dynamic,1> value;
76  // ------------------------------------
77  // 1) loop on elements
78  // ------------------------------------
79  result = Result(0);
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);
82  // ------------------------------------
83  // 1.1) locally compute values
84  // ------------------------------------
85  // locally evaluate int_K expr dx
86  // with loop on a quadrature formulae
87  // and accumulate in result
88  expr.evaluate (omega, K, value);
89  const Eigen::Matrix<T,Eigen::Dynamic,1>& w = pops.get_weight (omega, K);
90  check_macro (w.size() == value.size(),
91  "incompatible sizes w("<<w.size()<<") and value("<<value.size()<<")");
92  // TODO: DVT_BLAS1 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic> ri = value.transpose()*w;
93  // => compil pb with Eigen when Result is point, tensor, etc:
94  size_type loc_nnod = value.size();
95  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
96  result += value[loc_inod]*w[loc_inod];
97  }
98  }
99 #ifdef _RHEOLEF_HAVE_MPI
101  result = mpi::all_reduce (omega.geo_element_ownership(0).comm(), result, std::plus<Result>());
102  }
103 #endif // _RHEOLEF_HAVE_MPI
104 }
105 
106 }}// namespace rheolef::details
107 #endif // _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
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::integrate_option::get_order
size_t get_order() const
Definition: integrate_option.h:242
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::quadrature::set_order
void set_order(size_type order)
Definition: quadrature_rep.cc:284
rheolef::value
rheolef::std value
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::is_distributed
Definition: distributed.h:36
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::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef::quadrature::set_family
void set_family(family_type ft)
Definition: quadrature_rep.cc:292
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::details::field_expr_v2_value_assembly
void field_expr_v2_value_assembly(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result &result)
Definition: field_expr_value_assembly.h:55
size_type
field::size_type size_type
Definition: branch.cc:425
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