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
field::size_type size_type
Definition: branch.cc:425
see the Float page for the full documentation
size_type nnod(reference_element hat_K) const
characteristic_rep(const field_basic< T, M > &dh)
field_basic< T, M > _dh
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
see the disarray page for the full documentation
Definition: disarray.h:459
see the distributor page for the full documentation
Definition: distributor.h:62
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
size_type size(size_type iproc) const
Definition: distributor.h:163
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
const space_type & get_space() const
Definition: field.h:300
const geo_type & get_geo() const
Definition: field.h:301
void dis_dof_update() const
Definition: field.cc:410
see the geo_element page for the full documentation
Definition: geo_element.h:102
const quadrature< T > & get_quadrature() const
const basis_on_pointset< T > & get_basis_on_pointset() const
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
const quadrature_option & get_options() const
see the reference_element page for the full documentation
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
integrate_option quadrature_option
size_t size_type
Definition: basis_get.cc:76
distributed
Definition: asr.cc:228
point_basic< T >
Definition: piola_fem.h:135
#define warning_macro(message)
Definition: dis_macros.h:53
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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
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)
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)
This file is part of Rheolef.
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
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)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
disarray< point_basic< T >, M > _hat_y
disarray< point_basic< T >, M > _yq
point_basic< T > F
Definition: piola.h:79
Expr1::memory_type M
Definition: vec_expr_v2.h:416