Rheolef  7.1
an efficient C++ finite element environment
characteristic.cc
Go to the documentation of this file.
1 #include "rheolef/characteristic.h"
22 #include "rheolef/fem_on_pointset.h"
23 #include "rheolef/rheostream.h"
24 #include "rheolef/band.h"
25 #include "rheolef/field_evaluate.h"
26 
27 // TODO: replace piola_on_quad aka basis_on_pointset by piola_on_poinset
28 
29 // ===========================================================================
30 // allocator
31 // ===========================================================================
32 namespace rheolef {
33 
34 template<class T, class M>
36  : _dh(dh),
37  _coq_map()
38 {
39  // automatically add the "boundary" domain, if not already present:
40  _dh.get_geo().boundary();
41 }
42 
43 } // namespace rheolef
44 // ===========================================================================
45 // externs
46 // ===========================================================================
47 namespace rheolef { namespace details {
48 
49 template<class T, class M>
50 void
52  const geo_basic<T,M>& omega,
53  const disarray<point_basic<T>,M>& x,
54  const disarray<geo_element::size_type,M>& ix2dis_ie, // x has been already localized in K
55  disarray<index_set,M>& ie2dis_ix, // K -> list of ix
56  disarray<point_basic<T>,M>& hat_y);
57 
58 // ===========================================================================
59 // integrate(omega,compose (uh,X)) returns a field:
60 // It is done in two passes:
61 // - the first one is symbolic: localize the moved quadrature points yq=xq+X.dh(xq)
62 // this is done one time for all
63 // - the second pass is valued: compte uh(yq) during integration
64 // ===========================================================================
65 template <class T, class M>
66 void
68  // input :
69  const space_basic<T,M>& Xh,
70  const field_basic<T,M>& dh,
71  const piola_on_pointset<T>& pops,
72  // output :
73  disarray<index_set,M>& ie2dis_ix,
74  disarray<point_basic<T>,M>& hat_y,
76 {
77  // ----------------------------------------------------------------------
78  // 1) set the quadrature formulae
79  // ----------------------------------------------------------------------
80  typedef typename space_basic<T,M>::size_type size_type;
81  const geo_basic<T,M>& omega = Xh.get_geo();
82  check_macro (omega == dh.get_geo(), "characteristic: unsupported different meshes for flow ("
83  << dh.get_geo().name() << ") and convected field (" << omega.name() << ")");
84  if (omega.order() > 1) {
85  warning_macro ("characteristic: high-order curved elements not yet supported (treated as first order)");
86  }
87  fem_on_pointset<T> d_fops;
88  d_fops.initialize (dh.get_space().get_constitution().get_basis(), pops);
89  // ----------------------------------------------------------------------
90  // 2) size of the the disarray of all quadrature nodes
91  // ----------------------------------------------------------------------
92  // NOTE: since these nodes are moved (convected) by the flow associated to the characteristic
93  // they can change from partition and go to an element managed by another proc
94  // thus, in order to group comms, a whole disarray of quadrature point is allocated:
95  // In sequential mode (or with only one proc), this can do in an element by element way
96  // with less memory.
97  size_type dim = omega.dimension();
98  size_type map_dim = omega.map_dimension();
99  size_type sz = 0;
100  size_type dis_sz = 0;
101  const basis_on_pointset<T>& piola_on_quad = pops.get_basis_on_pointset();
104  distributor ige_ownership = omega.sizes().ownership_by_variant [variant];
105  if (ige_ownership.dis_size() == 0) continue;
106  reference_element hat_K (variant);
107  size_type nq = piola_on_quad.nnod (hat_K);
108  sz += nq*ige_ownership.size();
109  dis_sz += nq*ige_ownership.dis_size();
110  }
111  distributor xq_ownership (dis_sz, omega.comm(), sz);
112  // ----------------------------------------------------------------------
113  // 3) build the disarray of all quadrature nodes; xq
114  // ----------------------------------------------------------------------
115  disarray<point_basic<T>,M> xq (xq_ownership);
116  std::vector<size_type> dis_inod;
117  for (size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
118  const geo_element& K = omega[ie];
119  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
120  for (size_type q = 0, nq = piola.size(); q < nq; q++, ixq++) {
121  xq[ixq] = piola[q].F;
122  }
123  }
124  // ----------------------------------------------------------------------
125  // 4) interpolate dh on the xq nodes: dq(i)=dh(xq(i))
126  // ----------------------------------------------------------------------
127  dh.dis_dof_update();
128  disarray<point_basic<T>,M> dq (xq_ownership);
129  disarray<size_type,M> ix2dis_ie (xq_ownership);
130  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> d_value;
131  for (size_type ixq = 0, ie = 0, ne = omega.size(); ie < ne; ie++) {
132  const geo_element& K = omega[ie];
133  field_evaluate (dh, d_fops, omega, K, d_value);
134  for (size_type q = 0, nq = d_value.size(); q < nq; q++, ixq++) {
135  dq[ixq] = d_value [q];
136  ix2dis_ie[ixq] = ie; // will be a guest for locate(yq)
137  }
138  }
139  // ----------------------------------------------------------------------
140  // 5) build the disarray of all moved quadrature nodes yq(i) = xq(i)+dq(i)
141  // with the constraint to remain inside Omega
142  // ----------------------------------------------------------------------
143  yq.resize (xq_ownership);
144  omega.trace_move (xq, dq, ix2dis_ie, yq);
145  // ----------------------------------------------------------------------
146  // 6.a) symbolic interpolation-like pass: localize the yq nodes
147  // ----------------------------------------------------------------------
148  interpolate_pass1_symbolic (omega, yq, ix2dis_ie, ie2dis_ix, hat_y);
149 }
150 
151 } // namespace details
152 // ========================================================================
153 // symbolic pass 1 is stored one time for all in the characteristic class
154 // ========================================================================
155 template<class T, class M>
156 const characteristic_on_quadrature<T,M>&
158  const space_basic<T,M>& Xh,
159  const field_basic<T,M>& dh,
160  const piola_on_pointset<T>& pops) const
161 {
162  // get the quadrature used for integration
163  check_macro (pops.has_quadrature(), "unexpected point set, expect quadrature point set");
164  const quadrature<T>& quad = pops.get_quadrature();
165  const quadrature_option& qopt = quad.get_options();
166 
167  check_macro (qopt.get_family() == quadrature_option::gauss_lobatto,
168  "integrate characteristics HINT: use Gauss-Lobatto quadrature)");
169  check_macro (qopt.get_order() != std::numeric_limits<quadrature_option::size_type>::max(),
170  "integrate characteristics HINT: invalid unset order");
171 
172  // search if already used & memorized
173  std::string label = qopt.get_family_name() + "(" + itos(qopt.get_order()) + ")";
174  typename map_type::const_iterator iter = _coq_map.find (label);
175  if (iter != _coq_map.end()) {
176  return (*iter).second;
177  }
178  // build a new one & memorize it
179  // search if already used & memorized
182  coq_r._qopt = qopt;
183  coq_r._quad = quad;
184  coq_r._piola_on_quad = pops.get_basis_on_pointset();
185  details::integrate_pass1_symbolic (Xh, dh, pops,
186  coq_r._ie2dis_ix,
187  coq_r._hat_y,
188  coq_r._yq);
189 
190  // output :
191  std::pair <typename map_type::iterator,bool> iter_b = _coq_map.insert (std::make_pair(label,coq));
192  typename map_type::iterator iter2 = iter_b.first;
193  return (*iter2).second;
194 }
195 // ----------------------------------------------------------------------------
196 // instanciation in library
197 // ----------------------------------------------------------------------------
198 #define _RHEOLEF_instanciation(T,M) \
199 template class characteristic_rep<T,M>; \
200 
201 _RHEOLEF_instanciation(Float,sequential)
202 #ifdef _RHEOLEF_HAVE_MPI
204 #endif // _RHEOLEF_HAVE_MPI
205 
206 }// namespace rheolef
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
rheolef::point_basic
Definition: point.h:87
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::details::interpolate_pass1_symbolic
void interpolate_pass1_symbolic(const geo_basic< T, M > &omega, const disarray< point_basic< T >, M > &x, const disarray< geo_element::size_type, M > &ix2dis_ie, disarray< index_set, M > &ie2dis_ix, disarray< point_basic< T >, M > &hat_y)
Definition: interpolate.cc:44
rheolef::piola::F
point_basic< T > F
Definition: piola.h:79
rheolef::characteristic_on_quadrature_rep::_hat_y
disarray< point_basic< T >, M > _hat_y
Definition: characteristic.h:146
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::characteristic_on_quadrature_rep
Definition: characteristic.h:118
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::piola_on_pointset::get_quadrature
const quadrature< T > & get_quadrature() const
Definition: piola_on_pointset.h:224
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::details::integrate_pass1_symbolic
void integrate_pass1_symbolic(const space_basic< T, M > &Xh, const field_basic< T, M > &dh, const piola_on_pointset< T > &pops, disarray< index_set, M > &ie2dis_ix, disarray< point_basic< T >, M > &hat_y, disarray< point_basic< T >, M > &yq)
Definition: characteristic.cc:67
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
rheolef::quadrature::get_options
const quadrature_option & get_options() const
Definition: quadrature_rep.cc:264
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::characteristic_on_quadrature_rep::_piola_on_quad
basis_on_pointset< T > _piola_on_quad
Definition: characteristic.h:144
rheolef::characteristic_rep::characteristic_rep
characteristic_rep(const field_basic< T, M > &dh)
Definition: characteristic.cc:35
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::characteristic_on_quadrature_rep::_yq
disarray< point_basic< T >, M > _yq
Definition: characteristic.h:147
rheolef::fem_on_pointset::initialize
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
Definition: fem_on_pointset.h:130
rheolef::field_basic
Definition: field.h:235
rheolef::characteristic_rep::get_pre_computed
const characteristic_on_quadrature< T, M > & get_pre_computed(const space_basic< T, M > &Xh, const field_basic< T, M > &dh, const piola_on_pointset< T > &pops) const
Definition: characteristic.cc:157
rheolef::characteristic_on_quadrature_rep::_qopt
quadrature_option _qopt
Definition: characteristic.h:142
rheolef::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::distributor::dis_size
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
rheolef::fem_on_pointset
Definition: fem_on_pointset.h:117
rheolef::smart_pointer_base< characteristic_on_quadrature_rep< T, M >, details::constructor_copy< characteristic_on_quadrature_rep< T, M > > >::data
const characteristic_on_quadrature_rep< T, M > & data() const
Definition: smart_pointer.h:266
rheolef::characteristic_on_quadrature
Definition: characteristic.h:150
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
rheolef::characteristic_rep::_dh
field_basic< T, M > _dh
Definition: characteristic.h:187
rheolef::characteristic_on_quadrature_rep::_quad
quadrature< T > _quad
Definition: characteristic.h:143
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::basis_on_pointset
Definition: basis_on_pointset.h:192
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::basis_on_pointset::nnod
size_type nnod(reference_element hat_K) const
Definition: basis_on_pointset.h:263
rheolef::field_evaluate
void field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:39
rheolef::characteristic_on_quadrature_rep::_ie2dis_ix
disarray< index_set, M > _ie2dis_ix
Definition: characteristic.h:145
size_type
field::size_type size_type
Definition: branch.cc:425
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::distributed
distributed
Definition: asr.cc:228
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
rheolef::piola_on_pointset::has_quadrature
bool has_quadrature() const
Definition: piola_on_pointset.h:217
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::quadrature
Definition: quadrature.h:186
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::piola_on_pointset::get_basis_on_pointset
const basis_on_pointset< T > & get_basis_on_pointset() const
Definition: piola_on_pointset.h:210
rheolef::piola_on_pointset::get_piola
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
Definition: piola_on_pointset.h:239
rheolef::piola
Definition: piola.h:67