Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
space_constitution_assembly.cc
Go to the documentation of this file.
1
21//
22// helper for dof numbering with space product and component access
23// continuation of space_constitution.cc
24//
25// contents:
26// 8. dos_idof, for assembly
27//
28#include "rheolef/geo_domain.h"
29#include "rheolef/space.h"
30#include "rheolef/space_numbering.h"
31#include "rheolef/space_mult.h"
32#include "rheolef/piola_util.h"
33#include "rheolef/basis_get.h"
34#include <sstream>
35
36namespace rheolef {
37
38// ============================================================================
39// 8. dis_idof, for assembly
40// ============================================================================
41template <class T, class M>
44{
45 // lazy evaluated on the fly:
46 if ( _loc_ndof [hat_K.variant()] != std::numeric_limits<size_type>::max()) {
47 return _loc_ndof [hat_K.variant()];
48 }
49 if (! is_hierarchical()) {
50 size_type loc_comp_ndof = get_basis().ndof (hat_K);
51 size_type n_comp = (!_is_hier ? 1 : size());
52 _loc_ndof [hat_K.variant()] = n_comp*loc_comp_ndof;
53 return _loc_ndof [hat_K.variant()];
54 }
55 // mixed multi-component:
56 size_type loc_ndof = 0;
57 for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
58 const space_constitution<T,M>& constit = *iter;
59 loc_ndof += constit.loc_ndof(hat_K); // recursive call
60 }
61 _loc_ndof [hat_K.variant()] = loc_ndof;
62 return loc_ndof;
63}
64// DG:
65// integrate ("internal_sides", jump(u)*average(v));
66// => in that case, have double dofs on internal sides
67template <class T, class M>
68static
69bool is_dg_not_broken (
70 const space_constitution_rep<T,M>& constit,
71 const geo_basic<T,M>& omega_K,
72 const geo_element& K)
73{
74 bool is_dg = (!constit.get_basis().is_continuous()
75 && omega_K.map_dimension() + 1 == constit.get_geo().map_dimension()
76 && omega_K.map_dimension() == K.dimension());
77 // K is a side on omega_K and geo_space is_broken: two domains have the "is_broken" flag: "sides" and "internal_sides"
78 bool is_broken = constit.get_geo().is_broken(); // e.g. HDG lagrange multipliers
79 bool res = (is_dg && !is_broken);
80 return res;
81}
82// HDG:
83// space Mh(omega["sides"],"P1d"); trial lambda(Mh);
84// integrate (omega, on_local_sides(lambda));
85// => in that case, have to loop on sides
86template <class T, class M>
87static
88bool is_broken (
89 const space_constitution_rep<T,M>& constit,
90 const geo_basic<T,M>& omega_K)
91{
92 // only two domains have the "is_broken" flag: "sides" and "internal_sides"
93 bool res = constit.get_geo().is_broken() && (omega_K.map_dimension() == constit.get_geo().map_dimension() + 1);
94 return res;
95}
96// HDG-POST:
97// space Mhs(omega,"P1d[sides]"); trial lambda_s(Mhs);
98// integrate (omega, on_local_sides(lambda_s));
99// => in that case, have to loop on sides with bi-valued lambda_s
100template <class T, class M>
101static
102bool is_dg_restricted_to_interface (
103 const space_constitution_rep<T,M>& constit,
104 const geo_basic<T,M>& omega_K,
105 const geo_element& K)
106{
107 return constit.get_basis().option().is_trace_n()
108 && constit.get_geo().get_background_geo().map_dimension() == omega_K.map_dimension()
109 && omega_K.map_dimension() == K.dimension() + 1;
110}
111// assembly loop on geo_elements K on a domain omega_K
112// give its corresponding element K1 in space.geo
113template <class T, class M>
114const geo_element&
116 const geo_basic<T,M>& space_geo,
117 const geo_basic<T,M>& omega_K,
118 const geo_element& K_in)
119{
120 bool is_broken = space_geo.is_broken() && (omega_K.map_dimension() == space_geo.map_dimension() + 1);
121 bool use_dom2bgd = (omega_K.variant() == geo_abstract_base_rep<T>::geo_domain &&
122 space_geo.variant() != geo_abstract_base_rep<T>::geo_domain);
123 bool use_bgd2dom = (omega_K.variant() != geo_abstract_base_rep<T>::geo_domain &&
124 space_geo.variant() == geo_abstract_base_rep<T>::geo_domain);
125 if (!is_broken && use_bgd2dom) { // get from background mesh:
126 return space_geo.bgd2dom_geo_element (K_in);
127 } else if (!is_broken && use_dom2bgd) { // get on background mesh:
128 return omega_K.dom2bgd_geo_element (K_in);
129 } else { // usual case: use K directly
130 return K_in;
131 }
132}
133// assembly_loc_ndof
134// assembly_dis_idof : take into account DG on internal sides with doubled loc_ndof
135template <class T, class M>
138 const geo_basic<T,M>& omega_K,
139 const geo_element& K_in) const
140{
141 if (! is_hierarchical()) {
142trace_macro("loc_ndof: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<"...");
143 const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
144 size_type loc_ndof = 0;
145 if (is_dg_not_broken (*this, omega_K, K_in)) {
146trace_macro("loc_ndof: is_dg_not_broken...");
147 // dg: check for internal sides with doubled dofs
148 size_type map_d = K.dimension();
149 size_type dis_ie0 = K.master(0),
150 dis_ie1 = K.master(1);
151 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
152 "unexpected isolated side K_in="<<K_in.name()<<K_in.dis_ie() << " in omega_K=" << omega_K.name()
153 << " associated to K="<<K.name()<<K.dis_ie() << " in space_geo=" << get_geo().name());
154 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
155 // boundary side K
156 const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
157 loc_ndof = this->loc_ndof(L);
158 } else {
159 // internal side K
160 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
161 const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
162 loc_ndof = this->loc_ndof(L0) + this->loc_ndof(L1);
163 }
164 } else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
165trace_macro("loc_ndof: is_dg_restricted_to_interface...");
166 // HDG-POST multiplier: check for boundary sides dofs
167 // example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
168 size_type map_d = K.dimension();
169 size_type dis_ie0 = K.master(0),
170 dis_ie1 = K.master(1);
171 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
172 "unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
173 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
174 // boundary side K
175 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
177 L0.get_side_informations (K, sid0);
178 loc_ndof = get_basis().local_ndof_on_side (L0, sid0);
179 } else {
180 // internal side K
181 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
182 const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
183 side_information_type sid0, sid1;
184 L0.get_side_informations (K, sid0);
185 L1.get_side_informations (K, sid1);
186 loc_ndof = get_basis().local_ndof_on_side (L0, sid0)
187 + get_basis().local_ndof_on_side (L1, sid1);
188 // TODO: how to manage e.g. lambda_hs["interface"] = 1; it involves bi-valued side-based space (Pkd(S))^2
189 fatal_macro ("HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
190 }
191 } else if (is_broken (*this, omega_K)) {
192trace_macro("loc_ndof: is_broken...");
193 // HDG multiplier: check for boundary sides dofs
194 // example: omega_K="square" & get_geo="square[sides]"
195 size_type space_map_d = get_geo().map_dimension();
196 check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
197 reference_element hat_K = K;
198 for (size_type isid = 0, nsid = hat_K.n_side(); isid < nsid; ++isid) {
199 reference_element hat_S = hat_K.side(isid);
200 loc_ndof += get_basis().ndof(hat_S);
201 }
202 } else {
203trace_macro("loc_ndof: usual...");
204 // usual case:
205 loc_ndof = this->loc_ndof (K);
206 }
207trace_macro("loc_ndof: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<": result="<<loc_ndof<<" done");
208 return loc_ndof;
209 }
210 // hybrid multi-component:
211trace_macro("loc_ndof_hier: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<"...");
212 check_macro (_is_hier, "unexpected hierarchical space");
213 size_type loc_ndof = 0;
214 for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
215 const space_constitution<T,M>& constit = *iter;
216 loc_ndof += constit.assembly_loc_ndof(omega_K, K_in); // recursive call
217 }
218trace_macro("loc_ndof_hier: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<": result="<<loc_ndof<<" done");
219 return loc_ndof;
220}
221template <class T, class M>
222void
224 const geo_basic<T,M>& omega_K,
225 const geo_element& K_in,
226 typename std::vector<geo_element::size_type>::iterator& dis_idof_t,
227 const distributor& hier_ownership,
228 const std::vector<distributor>& start_by_flattened_component,
229 size_type& i_flat_comp) const
230{
231 size_type space_map_d = get_geo().map_dimension();
232 size_type map_d = K_in.dimension();
233 if (! is_hierarchical()) {
234 const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
235 size_type loc_flat_comp_ndof = 0;
236 if (is_dg_not_broken (*this, omega_K, K_in)) {
237 // DG case:
238 size_type dis_ie0 = K.master(0),
239 dis_ie1 = K.master(1);
240 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
241 "unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
242 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
243 // boundary side K with one adjacent element
244 const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
245 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L, dis_idof_t);
246 loc_flat_comp_ndof = get_basis().ndof (L);
247 } else {
248 // internal side K with two adjacent elements
249 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
250 const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
251 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_t);
252 size_type loc_flat_comp_ndof0 = get_basis().ndof (L0);
253 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, dis_idof_t + loc_flat_comp_ndof0);
254 loc_flat_comp_ndof = loc_flat_comp_ndof0 + get_basis().ndof (L1);
255 }
256 } else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
257 // HDG-POST multiplier: check for boundary sides dofs
258 // example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
259 size_type dis_ie0 = K.master(0),
260 dis_ie1 = K.master(1);
261 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
262 "unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
263 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
264 // boundary side K with one adjacent element
265 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
267 L0.get_side_informations (K, sid0);
268 std::vector<size_type> dis_idof_L0;
269 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_L0);
270 Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0;
271 get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
272 for (size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
273 size_type loc_idof = loc_idof_L0 [loc_sid_idof];
274 dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
275 }
276 loc_flat_comp_ndof = loc_idof_L0.size();
277 } else {
278 // internal side K with two adjacent elements
279 const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
280 const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
281 side_information_type sid0, sid1;
282 L0.get_side_informations (K, sid0);
283 L1.get_side_informations (K, sid1);
284 std::vector<size_type> dis_idof_L0, dis_idof_L1;
285 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_L0);
286 space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, dis_idof_L1);
287 Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0, loc_idof_L1;
288 get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
289 get_basis().local_idof_on_side (L1, sid1, loc_idof_L1);
290 check_macro (loc_idof_L0.size() == loc_idof_L1.size(), "unexpected bi-valued side: unmatch sizes");
291 for (size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
292 size_type loc_idof = loc_idof_L0 [loc_sid_idof];
293 dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
294 dis_idof_t [loc_sid_ndof+loc_sid_idof] = dis_idof_L1 [loc_idof];
295 }
296 loc_flat_comp_ndof = loc_idof_L0.size() + loc_idof_L1.size();
297 fatal_macro ("HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
298 }
299 } else if (is_broken (*this, omega_K)) {
300 // HDG multiplier: check for boundary sides dofs
301 // example: omega_K="square" & get_geo="square[sides]"
302 check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
303 for (size_type isid = 0, nsid = K_in.n_subgeo(space_map_d); isid < nsid; ++isid) {
304 size_type dis_isid_in = (K_in.dimension() == 1) ? K_in[isid] : (K_in.dimension() == 2) ? K_in.edge(isid) : K_in.face(isid);
305 const geo_element& S_in = omega_K.dis_get_geo_element (space_map_d, dis_isid_in);
306 bool have_bgd_dom = (get_geo().variant() == geo_abstract_base_rep<T>::geo_domain && get_geo().name() != omega_K.name());
307 const geo_element* S_ptr = (!have_bgd_dom) ? &S_in : &(get_geo().bgd2dom_geo_element(S_in));
308 const geo_element& S = *S_ptr;
309 space_numbering::dis_idof (get_basis(), get_geo().sizes(), S, dis_idof_t + loc_flat_comp_ndof);
310 loc_flat_comp_ndof += get_basis().ndof (S);
311 }
312 } else {
313 // non-broken cases: dg or usual
314 // usual: non-DG case:
315 space_numbering::dis_idof (get_basis(), get_geo().sizes(), K, dis_idof_t);
316 loc_flat_comp_ndof = get_basis().ndof (K);
317 }
318 // shift dis_idof depending on i_flat_comp
319 for (size_type loc_flat_comp_idof = 0; loc_flat_comp_idof < loc_flat_comp_ndof; ++loc_flat_comp_idof) {
320 size_type dis_flat_comp_idof = dis_idof_t [loc_flat_comp_idof];
321 distributor comp_ownership = ownership();
322 size_type iproc = comp_ownership.find_owner (dis_flat_comp_idof);
323 size_type first_dis_flat_comp_idof = comp_ownership.first_index(iproc);
324 size_type first_dis_idof = hier_ownership.first_index(iproc);
325 check_macro (dis_flat_comp_idof >= first_dis_flat_comp_idof, "unexpected dis_comp_idof");
326 size_type flat_comp_idof = dis_flat_comp_idof - first_dis_flat_comp_idof;
327 size_type start_flat_comp_idof = start_by_flattened_component [i_flat_comp].size(iproc);
328 size_type idof = start_flat_comp_idof + flat_comp_idof;
329 size_type dis_idof = first_dis_idof + idof;
330 dis_idof_t [loc_flat_comp_idof] = dis_idof;
331 }
332 dis_idof_t += loc_flat_comp_ndof;
333 i_flat_comp++;
334 return;
335 }
336 // hierarchical case
337 for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
338 const space_constitution<T,M>& comp_constit = operator[] (i_comp);
339 comp_constit.data()._assembly_dis_idof_recursive (omega_K, K_in, dis_idof_t, hier_ownership, start_by_flattened_component, i_flat_comp); // recursive call
340 }
341}
342template <class T, class M>
343void
345 const geo_basic<T,M>& omega_K,
346 const geo_element& K_in,
347 std::vector<geo_element::size_type>& dis_idof_t) const
348{
349 size_type i_flat_comp = 0;
350 size_type loc_ndof = assembly_loc_ndof (omega_K, K_in);
351 dis_idof_t.resize (loc_ndof);
352 typename std::vector<geo_element::size_type>::iterator first_dis_idof_t = dis_idof_t.begin();
353 _assembly_dis_idof_recursive (omega_K, K_in, first_dis_idof_t, ownership(), _start_by_flattened_component, i_flat_comp);
354}
355// ----------------------------------------------------------------------------
356// instanciation in library
357// ----------------------------------------------------------------------------
359
360#ifdef _RHEOLEF_HAVE_MPI
362#endif // _RHEOLEF_HAVE_MPI
363
364} // namespace rheolef
see the distributor page for the full documentation
Definition distributor.h:69
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
abstract base interface class
Definition geo.h:248
virtual std::string name() const =0
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
size_type dimension() const
size_type master(bool i) const
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
see the reference_element page for the full documentation
reference_element side(size_type loc_isid) const
variant_type variant() const
const basis_basic< T > & get_basis() const
void assembly_dis_idof(const geo_basic< T, M > &dom, const geo_element &bgd_K, std::vector< geo_element::size_type > &dis_idof) const
void _assembly_dis_idof_recursive(const geo_basic< T, M > &dom, const geo_element &bgd_K, typename std::vector< geo_element::size_type >::iterator &dis_idof_t, const distributor &hier_ownership, const std::vector< distributor > &start_by_flattened_component, size_type &i_flat_comp) const
size_type loc_ndof(const reference_element &hat_K) const
hierarchy_type::const_iterator const_iterator
hierarchy_type::size_type size_type
const geo_basic< T, M > & get_geo() const
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
size_type loc_ndof(const reference_element &hat_K) const
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
#define trace_macro(message)
Definition dis_macros.h:111
#define fatal_macro(message)
Definition dis_macros.h:33
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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.
const geo_element & assembly2space_geo_element(const geo_basic< T, M > &space_geo, const geo_basic< T, M > &omega_K, const geo_element &K_in)