Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
characteristic.cc
Go to the documentation of this file.
1
21#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// ===========================================================================
32namespace rheolef {
33
34template<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// ===========================================================================
47namespace rheolef { namespace details {
48
49template<class T, class M>
50void
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// ===========================================================================
65template <class T, class M>
66void
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,
76{
77 // ----------------------------------------------------------------------
78 // 1) set the quadrature formulae
79 // ----------------------------------------------------------------------
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();
103 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
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// ========================================================================
155template<class T, class M>
156const 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() + "(" + std::to_string(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;
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) \
199template class characteristic_rep<T,M>; \
200
202#ifdef _RHEOLEF_HAVE_MPI
203_RHEOLEF_instanciation(Float,distributed)
204#endif // _RHEOLEF_HAVE_MPI
205
206}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
size_type nnod(reference_element hat_K) const
characteristic_rep(const 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:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type dis_size() const
global and local sizes
size_type size(size_type iproc) const
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
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)
the finite element space
Definition space.h:382
#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)
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)
This file is part of Rheolef.
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)
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