Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
domain_indirect_seq.cc
Go to the documentation of this file.
1
21#include "rheolef/domain_indirect.h"
22#include "rheolef/geo.h"
23#include "rheolef/rheostream.h"
24
25namespace rheolef {
26
27// ----------------------------------------------------------------------------
28// @brief build a domain from a list of element indexes
29// ----------------------------------------------------------------------------
30template<class M>
31void
33 const std::string& name,
34 size_type map_dim,
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
53template<class M>
54void
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// ----------------------------------------------------------------------------
98void
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// ----------------------------------------------------------------------------
115template<class U>
116idiststream&
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}
162template<class U>
163void
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}
206template <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// ----------------------------------------------------------------------------
222template
225 idiststream& ips,
226 const geo_rep<Float,sequential>& omega,
227 std::vector<index_set> *ball);
228
229template
230void
232 const geo_rep<Float,sequential>& omega,
234 d_tmp,
235 std::vector<index_set>* ball);
236
237template
238domain_indirect_basic<sequential>::domain_indirect_basic (
240 d_tmp,
241 const geo_basic<Float, sequential>& omega,
242 std::vector<index_set>* ball);
243
244#define _RHEOLEF_instanciation(M) \
245template class domain_indirect_base_rep<M>; \
246
247_RHEOLEF_instanciation(sequential)
248#ifdef _RHEOLEF_HAVE_MPI
249_RHEOLEF_instanciation(distributed)
250#endif // _RHEOLEF_HAVE_MPI
251
252} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the disarray page for the full documentation
Definition disarray.h:497
rep::base::const_iterator const_iterator
Definition disarray.h:503
rep::base::iterator iterator
Definition disarray.h:502
see the distributor page for the full documentation
Definition distributor.h:69
static const size_type decide
Definition distributor.h:83
void build_from_list(const std::string &name, size_type map_dim, const communicator &comm, const std::vector< size_type > &ie_list)
void build_union(const domain_indirect_base_rep< M > &a, const domain_indirect_base_rep< M > &b)
geo_element_indirect::size_type size_type
the finite element boundary domain
abstract base interface class
Definition geo.h:248
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
geo_element_indirect::orientation_type orientation_type
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
size_type size() const
reference_element::size_type size_type
variant_type variant() const
geo_element_indirect::shift_type shift_type
size_type dis_ie() const
sequential mesh representation
Definition geo.h:778
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
std::istream & is()
Definition diststream.h:400
void inplace_intersection(const index_set &b)
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
see the smart_pointer page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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:44
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.
STL namespace.
Definition sphere.icc:25