Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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/*
27let:
28I = int_omega expr(x) dx
29
30The integrals are evaluated over each element K of the domain omega
31by using a quadrature formulae given by iopt
32
33expr(x) is a nonlinear expression
34
35The integrals over K are transformed on the reference element with
36the piola transformation:
37F : hat_K ---> K
38 hat_x |--> x = F(hat_x)
39
40int_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
44On each K, the computation needs a loop in q.
45The expresion is represented by a tree at compile-time.
46The 'expr' is represented by field_expr_terminal<f> : it evaluation gives :
47value1(q) = expr(F(hat_xq))
48
49This approch generilizes to an expression tree.
50*/
51namespace rheolef { namespace details {
52
53template <class T, class M, class Expr, class Result>
54void
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
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
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)
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.