Rheolef  7.2
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
23 #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& expr0,
58  const integrate_option& iopt,
59  Result& result)
60 {
61  Expr expr = expr0; // so could expr.initialze()
62  typedef typename geo_basic<T,M>::size_type size_type;
63  // ------------------------------------
64  // 0) initialize the loop
65  // ------------------------------------
66  quadrature<T> quad;
67  check_macro (iopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
68  "integrate: the quadrature formulae order may be chosen (HINT: use qopt.set_order(n))");
69  quad.set_order (iopt.get_order());
70  quad.set_family (iopt.get_family());
72  pops.initialize (omega.get_piola_basis(), quad, iopt);
73  expr.initialize (pops, iopt);
74  expr.template valued_check<Result>();
75  size_type map_d = omega.map_dimension();
76  Eigen::Matrix<Result,Eigen::Dynamic,1> value;
77  // ------------------------------------
78  // 1) loop on elements
79  // ------------------------------------
80  result = Result(0);
81  for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
82  const geo_element& K = omega.get_geo_element (map_d, ie);
83  // ------------------------------------
84  // 1.1) locally compute values
85  // ------------------------------------
86  // locally evaluate int_K expr dx
87  // with loop on a quadrature formulae
88  // and accumulate in result
89  expr.evaluate (omega, K, value);
90  const Eigen::Matrix<T,Eigen::Dynamic,1>& w = pops.get_weight (omega, K);
91  check_macro (w.size() == value.size(),
92  "incompatible sizes w("<<w.size()<<") and value("<<value.size()<<")");
93  // TODO: DVT_BLAS1 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic> ri = value.transpose()*w;
94  // => compil pb with Eigen when Result is point, tensor, etc:
95  size_type loc_nnod = value.size();
96  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
97  result += value[loc_inod]*w[loc_inod];
98  }
99  }
100 #ifdef _RHEOLEF_HAVE_MPI
102  result = mpi::all_reduce (omega.geo_element_ownership(0).comm(), result, std::plus<Result>());
103  }
104 #endif // _RHEOLEF_HAVE_MPI
105 }
106 
107 }}// namespace rheolef::details
108 #endif // _RHEO_FIELD_EXPR_VALUE_ASSEMBLY_H
field::size_type size_type
Definition: branch.cc:430
see the geo_element page for the full documentation
Definition: geo_element.h:102
see the integrate_option page for the full documentation
family_type get_family() const
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
const Eigen::Matrix< T, Eigen::Dynamic, 1 > & get_weight(const geo_basic< T, M > &omega, const geo_element &K) const
void set_family(family_type ft)
void set_order(size_type order)
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void field_expr_v2_value_assembly(const geo_basic< T, M > &omega, const Expr &expr0, const integrate_option &iopt, Result &result)
This file is part of Rheolef.