Rheolef  7.2
an efficient C++ finite element environment
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 
51 namespace rheolef { namespace details {
52 
53 template<class FieldWdof, class FieldLazy, class SetPlusOp>
54 typename 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;
69  using band_type = band_basic<scalar_type,memory_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
Definition: geo_element.h:102
typename scalar_traits< value_type >::type scalar_type
typename Expr1::memory_type memory_type
rheolef::std type
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)")
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)
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
This file is part of Rheolef.