Rheolef  7.1
an efficient C++ finite element environment
domain_indirect_seq.cc
Go to the documentation of this file.
1 #include "rheolef/domain_indirect.h"
22 #include "rheolef/geo.h"
23 #include "rheolef/rheostream.h"
24 
25 namespace rheolef {
26 
27 // ----------------------------------------------------------------------------
28 // @brief build a domain from a list of element indexes
29 // ----------------------------------------------------------------------------
30 template<class M>
31 void
33  const std::string& name,
35  const communicator& comm,
36  const std::vector<size_type>& ie_list)
37 {
38  _name = name;
39  _map_dim = map_dim;
40  distributor dom_ownership (distributor::decide, comm, ie_list.size());
43  for (std::vector<size_type>::const_iterator iter = ie_list.begin(),
44  last = ie_list.end();
45  iter != last; iter++, oige_iter++) {
46  (*oige_iter).set (1, *iter);
47  }
48 }
49 // ----------------------------------------------------------------------------
50 // @brief build the union of two domains
51 // ----------------------------------------------------------------------------
52 // c := a union b with c=*this
53 template<class M>
54 void
58 {
59  check_macro (a._map_dim == b._map_dim,
60  "union: domains "<<a._name<<" and " << b._name << " have incompatible dimensions");
61  _name = a._name + "+" + b._name;
62  _map_dim = a._map_dim;
63  // use a map to perform the union efficiently
64  std::map<size_type,geo_element_indirect> c_map;
66  iter = a.begin(), last = a.end(); iter != last; ++iter) {
67  size_type ie = (*iter).index();
68  c_map [ie] = *iter;
69  }
71  iter = b.begin(), last = b.end(); iter != last; ++iter) {
72  size_type ie = (*iter).index();
73  c_map [ie] = *iter; // not inserted if already present
74  }
75  distributor c_ownership (distributor::decide, a.comm(), c_map.size());
78  for (typename std::map<size_type,geo_element_indirect>::const_iterator
79  iter = c_map.begin(), last = c_map.end(); iter != last; ++iter, ++oige_iter) {
80  *oige_iter = (*iter).second;
81  }
82 }
83 // ----------------------------------------------------------------------------
85 // ----------------------------------------------------------------------------
88  using namespace std;
89  ops << "domain" << endl
90  << base::_name << endl
91  << "2 " << base::_map_dim << " " << size() << endl;
93  return ops;
94 }
95 // ----------------------------------------------------------------------------
97 // ----------------------------------------------------------------------------
98 void
100  const geo_element& S,
101  const std::vector<index_set>& ball,
102  index_set& contains_S)
103 {
105  contains_S.clear();
106  if (S.size() == 0) return;
107  contains_S = ball[S[0]];
108  for (size_type i = 1; i < S.size(); i++) {
109  contains_S.inplace_intersection (ball[S[i]]);
110  }
111 }
112 // ----------------------------------------------------------------------------
114 // ----------------------------------------------------------------------------
115 template<class U>
116 idiststream&
118  idiststream& ips,
119  const geo_rep<U,sequential>& omega,
120  std::vector<index_set> *ball)
121 {
122  std::istream& is = ips.is();
123  if (!scatch(is,"\ndomain")) {
124  is.setstate (std::ios::badbit);
125  return ips;
126  }
127  size_type version, n_side;
128  is >> base::_name >> version >> base::_map_dim >> n_side;
129  check_macro (version <= 2, "unexpected domain format version="<<version);
130  if (version == 2 || base::_map_dim == 0) {
133  // check that elements are sorted by increasing variant
134  {
135  size_type prev_variant = 0;
136  for (size_type ioige = 0, noige = size(); ioige < noige; ioige++) {
137  size_type ige = oige (ioige).index();
138  const geo_element& S = omega.get_geo_element (base::_map_dim, ige);
139  check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
140  prev_variant = S.variant();
141  }
142  }
143  return ips;
144  }
145  check_macro (version == 1, "unexpected domain format version="<<version);
146  check_macro (base::_map_dim <= omega.dimension(), "unexpected domain dimension="
147  <<base::_map_dim<<" > geometry dimension = " << omega.dimension());
148  //
149  // old version=1 format: searh for index of sides
150  // => this computation is slow, so format 2 is preferred
151  // when using large 3d meshes with large boundary domains
152  // or with 3d regions defined as domains
153  //
154  geo_element_auto<> elem_init;
156  disarray_t;
157  disarray_t d_tmp (n_side, elem_init);
158  d_tmp.get_values (ips);
159  build_from_data (omega, d_tmp, ball);
160  return ips;
161 }
162 template<class U>
163 void
165  const geo_rep<U,sequential>& omega,
167  d_tmp,
168  std::vector<index_set>* ball)
169 {
171  disarray_t;
172 
174  if (d_tmp.size() == 0) return;
175  base::_map_dim = d_tmp[0].dimension(); // get map_dim from d_tmp disarray elements
176 
177  if (ball [base::_map_dim].size() == 0) {
178  ball [base::_map_dim].resize (omega.n_vertex());
179  for (size_type ige = 0, nge = omega.geo_element_ownership(base::_map_dim).size(); ige < nge; ige++) {
180  const geo_element& E = omega.get_geo_element (base::_map_dim,ige);
181  for (size_type iloc = 0; iloc < E.size(); iloc++) {
182  ball [base::_map_dim] [E[iloc]] += ige;
183  }
184  }
185  }
186  disarray_t::const_iterator q = d_tmp.begin();
187  size_type prev_variant = 0;
188  for (iterator_ioige p = ioige_begin(), last = ioige_end(); p != last; ++p, ++q) {
189  const geo_element& S = *q;
190  check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
191  prev_variant = S.variant();
192  index_set contains_S;
193  build_set_that_contains_S (S, ball[base::_map_dim], contains_S);
194  check_macro (contains_S.size() >= 1, "domain element not in mesh: S.dis_ie=" << S.dis_ie());
195  check_macro (contains_S.size() <= 1, "problem with domain element: S.dis_ie=" << S.dis_ie());
196  // now, choose the only one mesh (sub-)element matching domain element S
197  size_type i_side = *(contains_S.begin());
198  const geo_element& S_orig = omega.get_geo_element (base::_map_dim, i_side);
201  check_macro (S_orig.get_orientation_and_shift (S, orient, shift),
202  "problem with domain element: S.dis_ie=" << S.dis_ie());
203  (*p).set (orient, i_side, shift);
204  }
205 }
206 template <class U>
209  d_tmp,
210  const geo_basic<U, sequential>& omega,
211  std::vector<index_set>* ball)
213 {
214  check_macro (omega.variant() == geo_abstract_base_rep<U>::geo,
215  "build a domain for a geo_domain: not yet");
216  const geo_rep<U,sequential>& omega_data = dynamic_cast<const geo_rep<U,sequential>&>(omega.data());
217  data().build_from_data (omega_data, d_tmp, ball);
218 }
219 // ----------------------------------------------------------------------------
220 // instanciation in library
221 // ----------------------------------------------------------------------------
222 template
223 idiststream&
225  idiststream& ips,
226  const geo_rep<Float,sequential>& omega,
227  std::vector<index_set> *ball);
228 
229 template
230 void
232  const geo_rep<Float,sequential>& omega,
234  d_tmp,
235  std::vector<index_set>* ball);
236 
237 template
240  d_tmp,
241  const geo_basic<Float, sequential>& omega,
242  std::vector<index_set>* ball);
243 
244 #define _RHEOLEF_instanciation(M) \
245 template class domain_indirect_base_rep<M>; \
246 
247 _RHEOLEF_instanciation(sequential)
248 #ifdef _RHEOLEF_HAVE_MPI
250 #endif // _RHEOLEF_HAVE_MPI
251 
252 } // namespace rheolef
rheolef::domain_indirect_base_rep::build_from_list
void build_from_list(const std::string &name, size_type map_dim, const communicator &comm, const std::vector< size_type > &ie_list)
Definition: domain_indirect_seq.cc:32
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
put
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
Definition: branch.cc:500
rheolef::geo_element::get_orientation_and_shift
bool get_orientation_and_shift(const geo_element &S, orientation_type &orient, shift_type &shift) const
return orientation and shift between *this element and S
Definition: geo_element.cc:114
rheolef::domain_indirect_base_rep::build_union
void build_union(const domain_indirect_base_rep< M > &a, const domain_indirect_base_rep< M > &b)
Definition: domain_indirect_seq.cc:55
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::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:126
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::sequential
Definition: distributed.h:28
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef-config.version
version
Definition: rheolef-config.in:126
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_rep
sequential mesh representation
Definition: geo.h:778
rheolef::geo_element::shift_type
geo_element_indirect::shift_type shift_type
Definition: geo_element.h:135
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::geo_element_auto
Definition: geo_element.h:376
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::domain_indirect_rep< sequential >::iterator_ioige
base::iterator_ioige iterator_ioige
Definition: domain_indirect.h:158
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::domain_indirect_base_rep::size_type
geo_element_indirect::size_type size_type
Definition: domain_indirect.h:77
p
Definition: sphere.icc:25
rheolef::geo_abstract_base_rep
abstract base interface class
Definition: geo.h:248
rheolef::smart_pointer
see the smart_pointer page for the full documentation
Definition: smart_pointer.h:351
rheolef::index_set::clear
void clear()
Definition: index_set_header.icc:152
rheolef::geo_element::size_type
reference_element::size_type size_type
Definition: geo_element.h:125
a
Definition: diffusion_isotropic.h:25
rheolef::domain_indirect_rep< sequential >::size_type
base::size_type size_type
Definition: domain_indirect.h:157
rheolef::scatch
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:52
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
size_type
field::size_type size_type
Definition: branch.cc:425
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::domain_indirect_base_rep
Definition: domain_indirect.h:71
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::domain_indirect_basic
the finite element boundary domain
Definition: domain_indirect.h:335
rheolef::std
Definition: vec_expr_v2.h:402
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
rheolef::domain_indirect_rep
Definition: domain_indirect.h:148
rheolef::build_set_that_contains_S
void build_set_that_contains_S(const geo_element &S, const std::vector< index_set > &ball, index_set &contains_S)
builds a set of elements that all contain S.
Definition: domain_indirect_seq.cc:99