rheolef  6.5
domain_indirect_seq.cc
Go to the documentation of this file.
1 #include "rheolef/domain_indirect.h"
2 #include "rheolef/geo.h"
3 #include "rheolef/rheostream.h"
4 
5 namespace rheolef {
6 
7 template<class M>
9  const std::string& name,
10  size_type map_dim,
11  const communicator& comm,
12  const std::vector<size_type>& ie_list)
13  : array<geo_element_indirect,M>(),
14  _name(),
15  _map_dim()
16 {
17  build_from_list (name, map_dim, comm, ie_list);
18 }
19 template<class M>
20 void
22  const std::string& name,
23  size_type map_dim,
24  const communicator& comm,
25  const std::vector<size_type>& ie_list)
26 {
27  _name = name;
28  _map_dim = map_dim;
29  distributor dom_ownership (distributor::decide, comm, ie_list.size());
30  array<geo_element_indirect,M>::resize (dom_ownership);
31  typename array<geo_element_indirect,M>::iterator oige_iter = array<geo_element_indirect,M>::begin();
32  for (std::vector<size_type>::const_iterator iter = ie_list.begin(),
33  last = ie_list.end();
34  iter != last; iter++, oige_iter++) {
35  (*oige_iter).set (1, *iter);
36  }
37 }
38 template<class M>
39 void
43 {
45  "union: domains "<<a._name<<" and " << b._name << " have incompatible dimensions");
46  _name = a._name + "+" + b._name;
47  _map_dim = a._map_dim;
48  std::map<size_type,geo_element_indirect> c_map;
49  for (typename array<geo_element_indirect,M>::const_iterator
50  iter = a.begin(), last = a.end(); iter != last; ++iter) {
51  size_type ie = (*iter).index();
52  c_map [ie] = *iter;
53  }
54  for (typename array<geo_element_indirect,M>::const_iterator
55  iter = b.begin(), last = b.end(); iter != last; ++iter) {
56  size_type ie = (*iter).index();
57  c_map [ie] = *iter; // not inserted if already present
58  }
59  distributor c_ownership (distributor::decide, a.comm(), c_map.size());
60  array<geo_element_indirect,M>::resize (c_ownership);
61  typename array<geo_element_indirect,M>::iterator oige_iter = array<geo_element_indirect,M>::begin();
62  for (typename std::map<size_type,geo_element_indirect>::const_iterator
63  iter = c_map.begin(), last = c_map.end(); iter != last; ++iter, ++oige_iter) {
64  *oige_iter = (*iter).second;
65  }
66 }
69  using namespace std;
70  ops << "domain" << endl
71  << base::_name << endl
72  << "2 " << base::_map_dim << " " << size() << endl;
73  array<geo_element_indirect,sequential>::put_values (ops);
74  return ops;
75 }
76 void
78  const geo_element& S,
79  const std::vector<index_set>& ball,
80  index_set& contains_S)
81 {
83  contains_S.clear();
84  if (S.size() == 0) return;
85  contains_S = ball[S[0]];
86  for (size_type i = 1; i < S.size(); i++) {
87  contains_S.inplace_intersection (ball[S[i]]);
88  }
89 }
90 template<class U>
93  idiststream& ips,
94  const geo_rep<U,sequential>& omega,
95  std::vector<index_set> *ball)
96 {
97  std::istream& is = ips.is();
98  if (!scatch(is,"\ndomain")) {
99  is.setstate (std::ios::badbit);
100  return ips;
101  }
102  size_type version, n_side;
103  is >> base::_name >> version >> base::_map_dim >> n_side;
104  check_macro (version <= 2, "unexpected domain format version="<<version);
105  if (version == 2 || base::_map_dim == 0) {
106  array<geo_element_indirect,sequential>::resize (n_side);
107  array<geo_element_indirect,sequential>::get_values (ips);
108  {
109  size_type prev_variant = 0;
110  for (size_type ioige = 0, noige = size(); ioige < noige; ioige++) {
111  size_type ige = oige (ioige).index();
112  const geo_element& S = omega.get_geo_element (base::_map_dim, ige);
113  check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
114  prev_variant = S.variant();
115  }
116  }
117  return ips;
118  }
119  check_macro (version == 1, "unexpected domain format version="<<version);
120  check_macro (base::_map_dim <= omega.dimension(), "unexpected domain dimension="
121  <<base::_map_dim<<" > geometry dimension = " << omega.dimension());
122  // => this computation is slow, so format 2 is prefered
124  geo_element_auto<heap_allocator<size_type> > elem_init (alloc);
125  typedef array<geo_element_auto<heap_allocator<size_type> >,sequential, heap_allocator<size_type> >
126  array_t;
127  array_t d_tmp (n_side, elem_init, alloc);
128  d_tmp.get_values (ips);
129  build_from_data (omega, d_tmp, ball);
130  return ips;
131 }
132 template<class U>
133 void
135  const geo_rep<U,sequential>& omega,
137  d_tmp,
138  std::vector<index_set>* ball)
139 {
140  typedef array<geo_element_auto<heap_allocator<size_type> >,sequential, heap_allocator<size_type> >
141  array_t;
142 
143  array<geo_element_indirect,sequential>::resize (d_tmp.size());
144  if (d_tmp.size() == 0) return;
145  base::_map_dim = d_tmp[0].dimension(); // get map_dim from d_tmp array elements
146 
147  if (ball [base::_map_dim].size() == 0) {
148  ball [base::_map_dim].resize (omega.n_vertex());
149  for (size_type ige = 0, nge = omega.geo_element_ownership(base::_map_dim).size(); ige < nge; ige++) {
150  const geo_element& E = omega.get_geo_element (base::_map_dim,ige);
151  for (size_type iloc = 0; iloc < E.size(); iloc++) {
152  ball [base::_map_dim] [E[iloc]] += ige;
153  }
154  }
155  }
156  array_t::const_iterator q = d_tmp.begin();
157  size_type prev_variant = 0;
158  for (iterator_ioige p = ioige_begin(), last = ioige_end(); p != last; ++p, ++q) {
159  const geo_element& S = *q;
160  check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
161  prev_variant = S.variant();
162  index_set contains_S;
163  build_set_that_contains_S (S, ball[base::_map_dim], contains_S);
164  check_macro (contains_S.size() >= 1, "domain element not in mesh: " << S);
165  check_macro (contains_S.size() <= 1, "problem with domain element: " << S);
166  size_type i_side = *(contains_S.begin());
167  const geo_element& S_orig = omega.get_geo_element (base::_map_dim, i_side);
170  check_macro (S_orig.get_orientation_and_shift (S, orient, shift),
171  "problem with domain element: " << S);
172  (*p).set (orient, i_side, shift);
173  }
174 }
175 template <class U>
178  d_tmp,
179  const geo_basic<U, sequential>& omega,
180  std::vector<index_set>* ball)
182 {
183  check_macro (omega.variant() == geo_abstract_base_rep<U>::geo,
184  "build a domain for a geo_domain: not yet");
185  const geo_rep<U,sequential>& omega_data = dynamic_cast<const geo_rep<U,sequential>&>(omega.data());
186  data().build_from_data (omega_data, d_tmp, ball);
187 }
188 template
191  idiststream& ips,
192  const geo_rep<Float,sequential>& omega,
193  std::vector<index_set> *ball);
194 
195 template
196 void
198  const geo_rep<Float,sequential>& omega,
200  d_tmp,
201  std::vector<index_set>* ball);
202 
203 template
206  d_tmp,
207  const geo_basic<Float, sequential>& omega,
208  std::vector<index_set>* ball);
209 
210 #define _RHEOLEF_instanciation(M) \
211 template class domain_indirect_base_rep<M>; \
212 
213 _RHEOLEF_instanciation(sequential)
214 #ifdef _RHEOLEF_HAVE_MPI
216 #endif // _RHEOLEF_HAVE_MPI
217 
218 } // namespace rheolef
219