Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_wdof_convert.h
Go to the documentation of this file.
1#ifndef _RHEO_FIELD_WDOF_CONVERT_H
2#define _RHEO_FIELD_WDOF_CONVERT_H
23// convert to field an un-assembled lazy_field expression
24// AUTHOR: Pierre.Saramito@imag.fr
25// DATE: 6 april 2020
26
27/*
28 Let:
29 uh = interpolate (Xh, expr);
30 or
31 l(v) = int_omega expr(v) dx
32 with v in Xh.
33
34 The expression is evaluated by a loop for each element K of the domain.
35
36 IMPLEMENTATION NOTE:
37 It writes
38 lh = integrate(omega, expr(v), iopt);
39 and the support for
40 lh += integrate(omega, expr(v), iopt);
41 lh -= integrate(omega, expr(v), iopt);
42 is provided with the "my_set_plus_op" argument, that is replaced
43 by details::generic_set_minus_op for -=
44
45 Supports for integration on a band or an usual domain.
46*/
47#include "rheolef/field_lazy.h"
48#include "rheolef/field_wdof.h"
49#include "rheolef/band.h"
50
51namespace rheolef { namespace details {
52
53template<class FieldWdof, class FieldLazy, class SetPlusOp>
54typename std::enable_if<
55 has_field_lazy_interface<FieldLazy>::value
56 && has_field_wdof_interface<FieldWdof>::value
57 ,FieldWdof&
58>::type
60 const FieldLazy& expr0,
61 const SetPlusOp& my_set_plus_op,
62 FieldWdof& uh)
63{
64 using size_type = typename FieldWdof::size_type;
65 using scalar_type = typename FieldWdof::scalar_type;
66 using memory_type = typename FieldWdof::memory_type;
67 using geo_type = typename FieldWdof::geo_type;
68 using space_type = typename FieldWdof::space_type;
70 FieldLazy expr = expr0; // so could call expr.initialize() which is non-const
71 // ------------------------------------
72 // 1) initialize the loop
73 // ------------------------------------
74 const space_type& Xh = expr.get_space();
75 check_macro (uh.get_space() == Xh, "invalid spaces");
76 geo_type omega = expr.get_geo();
77 expr.initialize (omega);
78 std::vector<size_type> dis_idx;
79 Eigen::Matrix<scalar_type,Eigen::Dynamic,1> lk;
80 size_type d = omega.dimension();
81 size_type map_d = omega.map_dimension();
82 if (Xh.get_constitution().is_discontinuous()) {
83 Xh.get_constitution().neighbour_guard(); // TODO: obsolete ?
84 }
85 bool is_on_band = expr.is_on_band();
86 band_type gh;
87 if (is_on_band) gh = expr.get_band();
88 // ------------------------------------
89 // 2) loop on elements
90 // ------------------------------------
91 for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
92 const geo_element& K = omega.get_geo_element (map_d, ie);
93 // ------------------------------------
94 // 2.1) locally compute dofs
95 // ------------------------------------
96 if (! is_on_band) {
97 Xh.get_constitution().assembly_dis_idof (omega, K, dis_idx);
98 } else { // on a band
99 size_type L_ie = gh.sid_ie2bnd_ie (ie);
100 const geo_element& L = gh.band() [L_ie];
101 Xh.dis_idof (L, dis_idx);
102 }
103 // ------------------------------------
104 // 2.2) locally compute values
105 // ------------------------------------
106 expr.evaluate (omega, K, lk);
107 // -----------------------------------------
108 // 2.3) assembly local lk in global field lh
109 // -----------------------------------------
110 check_macro (dis_idx.size() == size_type(lk.size()),
111 "incompatible sizes dis_idx("<<dis_idx.size()<<") and lk("<<lk.size()<<")");
112 for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
113 const scalar_type& value = lk [loc_idof];
114 size_type dis_idof = dis_idx [loc_idof];
115 my_set_plus_op (uh.dis_dof_entry (dis_idof), value);
116 }
117 }
118 // -----------------------------------------
119 // 3) finalize the assembly process
120 // -----------------------------------------
121 uh.dis_dof_update (generic_set_plus_op());
122 return uh;
123}
124
125} // namespace details
126
127}// namespace rheolef
128#endif // _RHEO_FIELD_WDOF_CONVERT_H
field::size_type size_type
Definition branch.cc:430
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the geo_element page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
std::enable_if< has_field_lazy_interface< FieldLazy >::value &&has_field_wdof_interface< FieldWdof >::value, FieldWdof & >::type convert_lazy2wdof(const FieldLazy &expr0, const SetPlusOp &my_set_plus_op, FieldWdof &uh)
This file is part of Rheolef.