Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_seq_get_bamg.cc
Go to the documentation of this file.
1
21//
22// geo: bamg input
23//
24// author: Pierre.Saramito@imag.fr
25//
26// 5 march 2012
27//
28#include "rheolef/geo.h"
29#include "rheolef/rheostream.h"
30#include "rheolef/iorheo.h"
31using namespace std;
32namespace rheolef {
33
34static
35void
36bamg_load_element (std::istream& is, size_t variant, geo_element_auto<>& K)
37{
38 K.reset (variant, 1);
39 for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
40 is >> K[iloc];
41 K[iloc]--;
42 }
43}
44template <class T>
45static
46void
47build_domains (
48 size_t dom_dim,
49 const disarray<geo_element_auto<>, sequential>& elt,
50 const disarray<size_t, sequential>& elt_bamg_dom_id,
51 const std::vector<std::string>& dom_name,
52 const geo_basic<T,sequential>& omega,
53 vector<index_set>* ball)
54{
55 using namespace std;
56 typedef typename geo_basic<T,sequential>::size_type size_type;
57 typedef pair<size_type,size_type> pair_t;
58 typedef geo_element_auto<> element_t;
59 typedef disarray<element_t, sequential> e_array_t;
60 typedef disarray<size_type, sequential> i_array_t;
61 typedef map<size_type, pair_t, less<size_type> > map_t;
62
63 // -------------------------------------------------------------------
64 // 1) reduce the bamg domain id to [0:ndom[ by counting domain id
65 // -------------------------------------------------------------------
66 // map[bamg_id] will gives idom and the number of its elements
67 map_t bamg_id2idom; // TODO (less<size_type>());
68 size_type idom = 0;
69 for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
70 size_type bamg_id = elt_bamg_dom_id [ie];
71 typename map_t::iterator iter = bamg_id2idom.find (bamg_id);
72 if (iter != bamg_id2idom.end()) {
73 // here is a previous bamg_dom_id: increment counter
74 ((*iter).second.second)++;
75 continue;
76 }
77 // here is a new bamg_dom_id: associated to idom and elt counter=1
78 bamg_id2idom.insert (pair<size_type,pair_t>(bamg_id, pair_t(idom,1)));
79 idom++;
80 }
81 size_type ndom = bamg_id2idom.size();
82 // -------------------------------------------------------------------
83 // 2) check that ndom matches the domain name disarray size
84 // -------------------------------------------------------------------
85 if (ndom != dom_name.size()) {
86 warning_macro ("geo::get_bamg: "
87 << dom_name.size() << " domain name(s) founded while "
88 << ndom << " bamg "<<dom_dim<<"d domain(s) are defined");
89 cerr << endl;
90 error_macro("HINT: check domain name file (.dmn)");
91 }
92 // -------------------------------------------------------------------
93 // 3) create domain disarray: loop on elements
94 // -------------------------------------------------------------------
95 element_t dummy_elt (reference_element::p, 0);
96 vector<e_array_t> domain (ndom, e_array_t(0, dummy_elt));
97 vector<size_type> counter (ndom, 0);
98 for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
99 size_type bamg_dom_id = elt_bamg_dom_id [ie];
100 size_type idom = bamg_id2idom [bamg_dom_id].first;
101 if (domain[idom].size() == 0) {
102 size_type size = bamg_id2idom [bamg_dom_id].second;
103 domain[idom] = e_array_t (size, dummy_elt);
104 }
105 size_type dom_ie = counter[idom];
106 domain[idom][dom_ie] = elt[ie];
107 counter[idom]++;
108 }
109 // -------------------------------------------------------------------
110 // 3) insert domains in the mesh: loop on domains
111 // -------------------------------------------------------------------
112 for (typename map_t::const_iterator
113 iter = bamg_id2idom.begin(),
114 last = bamg_id2idom.end(); iter != last; ++iter) {
115 size_type bamg_id = (*iter).first;
116 size_type idom = (*iter).second.first;
117 size_type size = (*iter).second.second;
118 check_macro (idom < dom_name.size(), "invalid idom="<<idom<<" for domain name");
119 string name = dom_name [idom];
120
121 domain_indirect_basic<sequential> dom (domain[idom], omega, ball);
122 dom.set_name (name);
123 dom.set_map_dimension (dom_dim);
124 omega.insert_domain_indirect (dom);
125 }
126 size_type ndom2 = omega.n_domain_indirect();
127}
128template <class T>
129static
130void
131build_vertex_domains (
132 const disarray<size_t, sequential>& edg_bdr_bamg_dom_id,
133 const disarray<size_t, sequential>& vert_bamg_dom_id,
134 const std::vector<std::string>& dom_name,
135 const geo_basic<T,sequential>& omega,
136 vector<index_set>* ball)
137{
138 if (dom_name.size() == 0) return;
139 typedef geo_element_auto<> element_t;
140 typedef disarray<element_t, sequential> v_array_t;
141 element_t dummy_elt (reference_element::p, 0);
142 // -------------------------------------------------------------------
143 // 1) list all vertex domain id
144 // -------------------------------------------------------------------
145 std::set<size_t> vert_id;
146 for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
147 size_t dom_id = vert_bamg_dom_id [iv];
148 if (dom_id == 0) continue;
149 vert_id.insert (dom_id);
150 }
151 // -------------------------------------------------------------------
152 // 2) omit vertex that are marked with an edge id
153 // -------------------------------------------------------------------
154 for (size_t iedg_bdr = 0, nedg_bdr = edg_bdr_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
155 size_t dom_id = edg_bdr_bamg_dom_id [iedg_bdr];
156 vert_id.erase (dom_id);
157 }
158 check_macro (vert_id.size() == dom_name.size(),
159 "unexpected VertexDomainNames with "<<dom_name.size()
160 <<" domain names while the mesh provides " << vert_id.size()
161 <<" vertex labels");
162 // -------------------------------------------------------------------
163 // 3) loop on vertex domain and insert it in mesh
164 // -------------------------------------------------------------------
165 size_t idom = 0;
166 for (std::set<size_t>::const_iterator iter = vert_id.begin(); iter != vert_id.end(); ++iter, ++idom) {
167 size_t id = *iter;
168 string name = dom_name[idom];
169 size_t dom_size = 0;
170 for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
171 if (vert_bamg_dom_id [iv] == id) dom_size++;
172 }
173 v_array_t vert_list (dom_size, dummy_elt);
174 for (size_t iv = 0, iv_dom = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
175 if (vert_bamg_dom_id [iv] != id) continue;
176 element_t& K = vert_list [iv_dom];
177 K.reset (reference_element::p, 1);
178 K[0] = iv;
179 iv_dom++;
180 }
181 domain_indirect_basic<sequential> dom (vert_list, omega, ball);
182 dom.set_name (name);
183 dom.set_map_dimension (0);
184 omega.insert_domain_indirect (dom);
185 }
186}
187template <class T>
188idiststream&
190{
191 using namespace std;
193 typedef typename geo_basic<T,sequential>::node_type node_type;
194
195 check_macro (ips.good(), "bad input stream for bamg");
196 istream& is = ips.is();
197 geo_header hdr;
198 hdr.dimension = 2;
199 hdr.map_dimension = 2;
200 hdr.order = 1;
201 hdr.dis_size_by_dimension [2] = 0;
202
203 // ------------------------------------------------------------------------
204 // 1) load data
205 // ------------------------------------------------------------------------
206 typedef geo_element_auto<> element_t;
207 typedef disarray<node_type, sequential> n_array_t;
208 typedef disarray<element_t, sequential> e_array_t;
209 typedef disarray<size_type, sequential> i_array_t;
210
211 n_array_t vertex;
212 i_array_t vert_bamg_dom_id;
213 i_array_t edg_bdr_bamg_dom_id;
214 i_array_t tri_bamg_dom_id;
215 i_array_t qua_bamg_dom_id;
216 e_array_t edg_bdr;
217 std::array<e_array_t, reference_element::max_variant> tmp_geo_element;
218 string label;
219 while (is.good() && label != "End") {
220 is >> label;
221 if (label == "Vertices") {
222 // ------------------------------------------------------------------------
223 // 1.1) load the coordinates
224 // Vertices <np>
225 // {xi yi dom_idx} i=0..np-1
226 // ------------------------------------------------------------------------
227 size_type nnod;
228 is >> nnod;
229 if (nnod == 0) {
230 warning_macro("empty bamg mesh file");
231 return ips;
232 }
233 hdr.dis_size_by_dimension [0] = nnod;
234 hdr.dis_size_by_variant [0] = nnod;
235 vertex.resize (nnod);
236 vert_bamg_dom_id = i_array_t (nnod, 0);
237 for (size_type inod = 0; inod < nnod; inod++) {
238 is >> vertex[inod][0] >> vertex[inod][1] >> vert_bamg_dom_id[inod];
239 }
240 check_macro (ips.good(), "bad input stream for bamg");
241 } else if (label == "Triangles") {
242 // ------------------------------------------------------------------------
243 // 2.1) load the connectivity
244 // Triangle <nt>
245 // {s0 s1 s2 dom_idx} j=0..nt-1
246 // ------------------------------------------------------------------------
247 size_type nt;
248 is >> nt;
249 hdr.dis_size_by_dimension [2] += nt;
251 element_t init_tri (reference_element::t, hdr.order);
252 tmp_geo_element [reference_element::t] = e_array_t (nt, init_tri);
253 tri_bamg_dom_id = i_array_t (nt, 0);
254 for (size_type it = 0; it < nt; it++) {
255 bamg_load_element (is, reference_element::t,
256 tmp_geo_element[reference_element::t][it]);
257 is >> tri_bamg_dom_id[it];
258 }
259 } else if (label == "Quadrilaterals") {
260 // ------------------------------------------------------------------------
261 // Quadrilaterals <nq>
262 // {s0 s1 s2 s3 dom_idx} j=0..nq-1
263 // ------------------------------------------------------------------------
264 size_type nq;
265 is >> nq;
266 hdr.dis_size_by_dimension [2] += nq;
268 element_t init_qua (reference_element::q, hdr.order);
269 tmp_geo_element [reference_element::q] = e_array_t (nq, init_qua);
270 qua_bamg_dom_id = i_array_t (nq, 0);
271 for (size_type iq = 0; iq < nq; iq++) {
272 bamg_load_element (is, reference_element::q,
273 tmp_geo_element[reference_element::q][iq]);
274 is >> qua_bamg_dom_id[iq];
275 }
276 } else if (label == "Edges") {
277 // ------------------------------------------------------------------------
278 // 2.3) load the boundary domains
279 // Edges <nedg>
280 // {s0 s1 dom_idx} j=0..nedg-1
281 // ------------------------------------------------------------------------
282 size_type nedg;
283 is >> nedg;
284 element_t init_edg (reference_element::e, hdr.order);
285 edg_bdr = e_array_t (nedg, init_edg);
286 edg_bdr_bamg_dom_id = i_array_t (nedg, 0);
287 for (size_type iedg = 0; iedg < nedg; iedg++) {
288 element_t& K = edg_bdr[iedg];
289 K.reset (reference_element::e, hdr.order);
290 for (size_type iloc = 0, nloc = 2; iloc < nloc; iloc++) {
291 is >> K[iloc];
292 K[iloc]--;
293 }
294 is >> edg_bdr_bamg_dom_id[iedg];
295 }
296 }
297 }
298 // ---------------------------------------------------------------
299 // 2) check rheolef extension to optional domain names
300 // ---------------------------------------------------------------
301 vector<string> vertice_domain_name;
302 vector<string> edg_dom_name;
303 vector<string> region_domain_name;
304 char c;
305 is >> ws >> c; // skip white and grab a char
306 // have "EdgeDomainNames" or "VertexDomainNames" ?
307 // bamg mesh may be followed by field data and such, so be carrefull...
308 while (c == 'E' || c == 'V' || c == 'R') {
309 is.unget(); // put char back
310 if (c == 'R') {
311 if (!scatch(is,"RegionDomainNames")) break;
312 size_type n_dom_region;
313 is >> n_dom_region;
314 region_domain_name.resize (n_dom_region);
315 for (size_type k = 0; k < n_dom_region; k++) {
316 is >> region_domain_name[k];
317 }
318 } else if (c == 'E') {
319 if (!scatch(is,"EdgeDomainNames")) break;
320 // ---------------------------------------------------------------
321 // get edge domain names
322 // ---------------------------------------------------------------
323 size_type n_dom_edge;
324 is >> n_dom_edge;
325 edg_dom_name.resize (n_dom_edge);
326 for (size_type k = 0; k < n_dom_edge; k++) {
327 is >> edg_dom_name[k];
328 }
329 } else if (c == 'V') {
330 if (!scatch(is,"VertexDomainNames")) break;
331 // ---------------------------------------------------------------
332 // get vertice domain names
333 // ---------------------------------------------------------------
334 size_type n_dom_vertice;
335 is >> n_dom_vertice;
336 vertice_domain_name.resize (n_dom_vertice);
337 for (size_type k = 0; k < n_dom_vertice; k++) {
338 is >> vertice_domain_name[k];
339 }
340 }
341 is >> ws >> c; // skip white and grab a char
342 }
343 // ------------------------------------------------------------------------
344 // 3) build & upgrade
345 // ------------------------------------------------------------------------
346 bool do_upgrade = true;
347 omega.build_from_data (hdr, vertex, tmp_geo_element, do_upgrade);
348
349 // ------------------------------------------------------------------------
350 // 4) get domain, until end-of-file
351 // ------------------------------------------------------------------------
352 vector<index_set> ball[4];
353 build_domains (1, edg_bdr, edg_bdr_bamg_dom_id, edg_dom_name, omega, ball);
354 build_vertex_domains (edg_bdr_bamg_dom_id, vert_bamg_dom_id, vertice_domain_name, omega, ball);
355 //TODO: region(d=2) domains
356 return ips;
357}
358// ----------------------------------------------------------------------------
359// instanciation in library
360// ----------------------------------------------------------------------------
362
363}// namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the disarray page for the full documentation
Definition disarray.h:497
generic mesh with rerefence counting
Definition geo.h:1089
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
std::istream & is()
Definition diststream.h:400
static const variant_type q
static const variant_type e
static const variant_type p
static const variant_type t
#define error_macro(message)
Definition dis_macros.h:49
#define warning_macro(message)
Definition dis_macros.h:53
const point vertex[n_vertex]
Definition edge.icc:67
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.
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
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
template idiststream & geo_get_bamg< Float >(idiststream &, geo_basic< Float, sequential > &)
STL namespace.
size_type map_dimension
Definition geo_header.h:40
size_type dis_size_by_dimension[4]
Definition geo_header.h:44
size_type dimension
Definition geo_header.h:39
size_type dis_size_by_variant[reference_element::max_variant]
Definition geo_header.h:43