2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
6
/// Rheolef is free software; you can redistribute it and/or modify
7
/// it under the terms of the GNU General Public License as published by
8
/// the Free Software Foundation; either version 2 of the License, or
9
/// (at your option) any later version.
11
/// Rheolef is distributed in the hope that it will be useful,
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14
/// GNU General Public License for more details.
16
/// You should have received a copy of the GNU General Public License
17
/// along with Rheolef; if not, write to the Free Software
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
/// =========================================================================
25
// Nicolas.Roquet@imag.fr
26
// Pierre.Saramito@imag.fr
30
// added label_interface and build_subregion for updating meshes for
31
// multi-region problems. J.Etienne@damtp.cam.ac.uk, 19/9/2006
33
// ============================================================================
35
// ============================================================================
37
#include "rheolef/georep.h"
38
#include "rheolef/iorheo.h"
39
#include "rheolef/rheostream.h" // i/o utility
40
#include "rheolef/tiny_element.h"
41
#include "rheolef/tensor.h"
42
#include "geo-connectivity.h"
46
// ============================================================================
48
// ============================================================================
53
typedef vector<int>::size_type size_type;
55
T& operator[](size_type i) { return _t[i]; }
56
const T& operator[](size_type i) const { return _t[i]; }
61
class bidimmensional_array {
63
typedef vector<int>::size_type size_type;
64
bidimmensional_array(size_type ni=0, size_type nj=0)
65
: _nj(nj), _t(ni*nj) {}
66
void resize (size_type ni, size_type nj) { _nj = nj; _t.resize(ni*nj); }
67
T* operator[](size_type i) { return _t.begin() + i*_nj; }
68
T* operator()(size_type i) { return operator[](i); }
69
T& operator()(size_type i, size_type j) { return _t[i*_nj+j]; }
74
// ============================================================================
75
// coordinate system mini-handler
76
// ============================================================================
78
georep::coordinate_system () const
80
return fem_helper::coordinate_system_name (_coord_sys);
83
georep::set_coordinate_system (const string& name)
85
_coord_sys = fem_helper::coordinate_system(name);
87
// ============================================================================
89
// ============================================================================
91
georep::georep(const georep& g)
93
Vector<geo_element>(g),
96
_serial_number(g._serial_number),
98
_coord_sys(g._coord_sys),
100
_domlist(g._domlist),
103
_save_double_precision(g._save_double_precision),
104
_version(g._version),
108
_localizer_initialized(false),
109
_boundary(g._boundary),
110
_boundary_initialized(g._boundary_initialized)
112
for (size_type i = 0; i < 4; i++)
113
_count_geo[i] = g._count_geo[i];
114
for (size_type i = 0; i < reference_element::max_size; i++)
115
_count_element[i] = g._count_element[i];
117
// ----------------------------------------------------------------------------
118
// cstor from file name
119
// ----------------------------------------------------------------------------
122
georep::reset_counters ()
124
fill(_count_geo, _count_geo+4, 0);
125
fill(_count_element, _count_element+reference_element::max_size, 0);
127
georep::georep(const string& a_file_name, const string& coord_sys_name)
129
Vector<geo_element>(),
131
_number(numeric_limits<size_type>::max()),
132
_serial_number(numeric_limits<size_type>::max()),
134
_coord_sys(fem_helper::coordinate_system(coord_sys_name)),
139
_save_double_precision(false),
144
_localizer_initialized(false),
146
_boundary_initialized()
149
string file_name = a_file_name;
150
irheostream is(file_name, "geo");
151
if (!is.good()||(file_name.find("field")!=string::npos))
153
warning_macro("\"" << file_name << "[.geo[.gz]]\" not found, trying for \"" << file_name << "[.[m]field[.gz]]\".");
154
is.open(file_name, "mfield");
155
if (!is.good()) is.open(file_name, "field");
156
check_macro(is.good(), "\"" << file_name << "[.geo[.gz]]\" not found.");
157
scatch(is,"\nfield");
160
is >> version >> ndof >> file_name;
161
is.open(file_name,"geo");
162
check_macro(is.good(), "\"" << file_name << "[.geo[.gz]]\" not found.");
164
string root_name = delete_suffix (delete_suffix(file_name, "gz"), "geo");
165
string base_name = get_basename (root_name);
166
string::size_type i = base_name.find_last_of ('+');
167
string::size_type j = base_name.find_last_of ('-');
168
string::size_type l = base_name.length();
172
string digits ="0123456789";
173
string tmp =string(base_name, i+1, j-i-1);
174
if (tmp.find_first_not_of(digits)==string::npos)
175
_number =atoi(tmp.c_str());
176
else // it is not a valid version number
178
warning_macro(tmp << " is not a valid mesh number in filename " << base_name);
183
tmp =string(base_name, j+1, l-j-1);
184
if (tmp.find_first_not_of(digits)==string::npos)
185
_serial_number =atoi(tmp.c_str());
188
warning_macro(tmp << " is not a valid mesh serial number in filename " << base_name);
189
i=l; j=l; _number=numeric_limits<size_type>::max();
192
base_name = string(base_name, 0, i);
194
is >> setbasename(base_name) >> *this; // load data
195
fem_helper::check_coord_sys_and_dimension (_coord_sys, _dim);
198
// ----------------------------------------------------------------------------
199
// Cstructor from point list and former geo
200
// ----------------------------------------------------------------------------
202
georep::georep (const nodelist_type& p, const georep& g)
204
Vector<geo_element>(g),
206
_number(numeric_limits<size_type>::max()),
207
_serial_number(numeric_limits<size_type>::max()),
209
_coord_sys(g._coord_sys),
211
_domlist(g._domlist),
214
_save_double_precision(g._save_double_precision),
215
_version(g._version),
219
_localizer_initialized(false),
220
_boundary(g._boundary),
221
_boundary_initialized(g._boundary_initialized)
223
for (size_type j = 0; j < _dim; j++) {
224
_xmin[j] = numeric_limits<Float>::max();
225
_xmax[j] = -numeric_limits<Float>::max();
227
for (size_type j = _dim+1; j < 3; j++) {
228
_xmin [j] = _xmax [j] = 0;
230
nodelist_type::const_iterator iter_p = p.begin();
231
nodelist_type::const_iterator last_p = p.end();
233
while (iter_p != last_p) {
234
const point& p = *iter_p++;
235
for (size_type j = 0 ; j < _dim; j++) {
236
_xmin[j] = min(p[j], _xmin[j]);
237
_xmax[j] = max(p[j], _xmax[j]);
241
if (g._number == numeric_limits<size_type>::max())
244
_number = g._number+1;
245
_serial_number = numeric_limits<size_type>::max();
247
for (size_type i = 0; i < 4; i++)
248
_count_geo[i] = g._count_geo[i];
249
for (size_type i = 0; i < reference_element::max_size; i++)
250
_count_element[i] = g._count_element[i];
252
// ----------------------------------------------------------------------------
254
// ----------------------------------------------------------------------------
256
georep::georep(const georep& g, const domain& d)
258
Vector<geo_element>(d.size()),
260
_number(numeric_limits<size_type>::max()),
261
_serial_number(numeric_limits<size_type>::max()),
263
_coord_sys(g._coord_sys),
268
_save_double_precision(g._save_double_precision),
273
_localizer_initialized(false),
275
_boundary_initialized()
277
g.may_have_version_2();
279
// step 0 : rename the mesh
281
_name = d.name() + "_from_" + g.name() ;
283
// step 1 : numbering of vertices
285
vector<bool> node_on_d (g.n_node(), false);
286
domain::const_iterator last_side = d.end();
287
for (domain::const_iterator i = d.begin(); i != last_side; i++) {
289
const geo_element& S = *i;
290
for (size_type j = 0; j < S.size(); j++)
291
node_on_d [S[j]] = true ;
293
size_type nvertex = 0;
294
vector<size_type> new_vertex_idx(g.n_node(), numeric_limits<size_type>::max());
295
for (size_type i = 0; i < g.n_node(); i++) {
297
new_vertex_idx[i] = nvertex++ ;
300
for (size_type j = 0; j < _dim; j++) {
301
_xmin[j] = numeric_limits<Float>::max();
302
_xmax[j] = -numeric_limits<Float>::max();
304
for (size_type j = _dim+1; j < 3; j++) {
305
_xmin [j] = _xmax [j] = 0;
308
size_type i_g_node = 0;
309
iterator_node iter_node = begin_node() ;
310
for (const_iterator_node iter_g_node = g.begin_node(); iter_g_node != g.end_node(); iter_g_node++, i_g_node++) {
311
if (node_on_d[i_g_node]) {
312
*iter_node = *iter_g_node;
313
point& p = *iter_node;
314
for (size_type j = 0 ; j < _dim; j++) {
315
_xmin[j] = min(p[j], _xmin[j]);
316
_xmax[j] = max(p[j], _xmax[j]);
321
// step 2 : copy geo_element from "d" into geo_element of *this
323
_count_geo [0] = n_node();
324
_count_element [0] = n_node();
327
iterator iter_elt = begin();
328
for (domain::const_iterator iter_side = d.begin();
329
iter_side != last_side;
330
iter_side++, K_idx++, iter_elt++) {
332
const geo_element& S = *iter_side;
333
geo_element& K = *iter_elt;
334
K.set_type(S.type()) ;
336
for (size_type i = 0; i < K.size(); i++) {
337
K[i] = new_vertex_idx [S[i]];
340
_count_geo [K.dimension()]++;
341
_count_element[K.type()]++;
343
build_connectivity (*this);
346
// ============================================================================
348
// ============================================================================
350
georep::build_subregion(const domain& start_from, const domain& dont_cross,
351
std::string name, std::string name_of_complement)
353
if (has_domain(name))
355
if (name_of_complement=="" || has_domain(name_of_complement))
357
warning_macro("Domain `" << name << "' already exists! Not overwriting.");
361
error_macro("Domain `" << name << "' already exists, but not `"<<name_of_complement<<"'!");
363
else if (name_of_complement!="" && has_domain(name_of_complement))
364
error_macro("Domain `" << name_of_complement << "' already exists, but not `"<<name<<"'!");
365
bool update_complement=(name_of_complement!="");
366
vector<set<size_type> > ball (n_node());
367
build_point_to_element_sets (begin(), end(), ball.begin());
368
// remember which are the vertices we don't want to cross
369
set< size_type > no_way_vertices;
370
domain::const_iterator iE=dont_cross.begin();
371
domain::const_iterator eE=dont_cross.end();
372
for ( ; iE!=eE; iE++)
373
for (size_type i=0; i<(*iE).size(); i++)
374
no_way_vertices.insert((*iE)[i]);
375
// a list to remember the furthest we went in our exploration so far
376
slist< size_type > frontier_dofs; // contains dofs
377
// ...and a set to mark elemnts we've seen already
378
set< size_type > elementa_cognita; // contains element indices
379
// take a vertex in start_from to start with:
380
domain::const_iterator i_start=start_from.begin();
381
domain::const_iterator i_start_end=start_from.end();
382
size_type v=(*i_start)[0];
383
while (no_way_vertices.count(v)!=0)
386
check_macro(i_start!=i_start_end, "No domain found between " << start_from.name() << " and "
387
<< dont_cross.name() );
390
frontier_dofs.push_front(v);
393
v=*(frontier_dofs.begin());
394
frontier_dofs.pop_front();
395
set<size_type>::const_iterator iK=ball[v].begin();
396
set<size_type>::const_iterator eK=ball[v].end();
398
if (elementa_cognita.count(*iK)==0)
400
elementa_cognita.insert(*iK);
401
for (size_type i=0; i<begin()[*iK].size(); i++)
402
if ( (begin()[*iK][i]!=v) && (no_way_vertices.count(begin()[*iK][i])==0) )
403
frontier_dofs.push_front(begin()[*iK][i]);
405
} while (!frontier_dofs.empty());
406
// region has been explored, now fill in domain structures for it and its complement
407
domain region, c_region;
408
region.set_name(name);
409
region.set_dimension(_dim);
410
c_region.set_name(name_of_complement);
411
c_region.set_dimension(_dim);
412
const_iterator iK=begin();
413
const_iterator eK=end();
414
for ( ; iK!=eK; iK++)
415
if (elementa_cognita.count((*iK).index())==0)
416
c_region.push_back(*iK);
418
region.push_back(*iK);
419
// insert regions in domain list
420
_domlist.push_back(region);
421
if (update_complement) _domlist.push_back(c_region);
424
// ============================================================================
425
// hazel initialization
426
// ============================================================================
428
georep::init_localizer (const domain& boundary, Float tolerance, int list_size) const
430
if (_localizer_initialized)
432
switch (dimension()) {
434
_h1d.initialize (*this);
437
_h.grow (begin(), end(), _x.begin(), _x.end(),
438
boundary.begin(), boundary.end(),
439
size(), dimension(), xmin(), xmax(),
440
tolerance, list_size);
443
_tree.initialize (map_dimension(), xmin(), xmax(),
444
n_vertex(), begin_node(), size(), begin());
446
default : error_macro("localizer not defined in " << dimension() << "D");
448
_localizer_initialized = true;
450
// ============================================================================
452
// ============================================================================
458
georep::has_domain (const string& domname) const
460
domlist_type::const_iterator iter = _domlist.begin();
461
domlist_type::const_iterator last = _domlist.end();
462
while (iter != last) {
463
const domain& dom = *iter;
464
if (domname == dom.name())
471
georep::operator[] (const string& domname) const
473
domlist_type::const_iterator iter = _domlist.begin();
474
domlist_type::const_iterator last = _domlist.end();
475
while (iter != last) {
476
const domain& dom = *iter;
477
if (domname == dom.name())
481
error_macro ("geo: undefined domain `" << domname << "'");
483
return *iter;// for no warning at compile time
485
// ============================================================================
487
// ============================================================================
496
build_connectivity (*this);
499
// ============================================================================
501
// ============================================================================
504
georep::map_dimension() const
506
for (size_type d = 3; d >= 0; d--) {
507
if (_count_geo [d] != 0) {
511
// has no elements...
518
may_have_version_2();
519
const_iterator_node p = begin_node();
520
const_iterator iter_elt = begin();
521
const_iterator last_elt = end();
522
Float hm2 = numeric_limits<Float>::max();
523
while (iter_elt != last_elt) {
524
const geo_element& K = *iter_elt++;
525
for (size_type e = 0; e < K.n_edge(); e++) {
526
size_type i1 = K.subgeo_vertex (1, e, 0);
527
size_type i2 = K.subgeo_vertex (1, e, 1);
528
hm2 = min(hm2, norm2(p[i1] - p[i2]));
536
may_have_version_2();
537
const_iterator_node p = begin_node();
538
const_iterator iter_elt = begin();
539
const_iterator last_elt = end();
541
while (iter_elt != last_elt) {
542
const geo_element& K = *iter_elt++;
543
for (size_type e = 0; e < K.n_edge(); e++) {
544
size_type i1 = K.subgeo_vertex (1, e, 0);
545
size_type i2 = K.subgeo_vertex (1, e, 1);
546
hm2 = max(hm2, norm2(p[i1] - p[i2]));
551
// =================================================================
552
// hatter and dehatter:
553
// F : hat(K) --> K is a nonlinear transformation
554
// that depends upon the degree of the polynomial basis
555
// => depend on the "space" class
556
// => cannot be handled by the "geo" class
557
// => hatter and dehatter may move to the "piola" class (TODO)
558
// =================================================================
560
georep::hatter (const point& x, size_type K_idx) const
562
const geo_element& K = element(K_idx);
563
switch(begin()[K_idx].type()) {
564
case geo_element::e: {
565
double a = begin_node()[begin()[K_idx][0]][0];
566
double b = begin_node()[begin()[K_idx][1]][0];
567
return meshpoint(K_idx, point((x[0] - a)/(b-a)));
569
case geo_element::t: {
570
point a = begin_node()[begin()[K_idx][0]];
571
point b = begin_node()[begin()[K_idx][1]];
572
point c = begin_node()[begin()[K_idx][2]];
574
Float t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
575
Float t11 = -a[0]+xx[0];
576
Float t15 = -a[1]+xx[1];
578
return meshpoint(K_idx,
579
point((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
580
(b[1]-a[1])*t9*t11-(b[0]-a[0])*t9*t15,
584
case geo_element::T: {
585
const point& a = vertex(K[0]);
586
const point& b = vertex(K[1]);
587
const point& c = vertex(K[2]);
588
const point& d = vertex(K[3]);
591
for (size_t i = 0; i < 3; i++) {
598
bool is_singular = ! invert_3x3 (A, inv_A);
599
check_macro(!is_singular, "hatter: singular transformation in element: " << K.index());
600
point hat_x = inv_A*ax;
601
return meshpoint(K_idx, hat_x);
603
default : error_macro("hatter: element not yet implemented: " << K.name());
607
georep::dehatter (const meshpoint& S, bool is_a_vector) const
609
switch(begin()[S.element()].type())
611
case geo_element::e: {
612
double a = begin_node()[begin()[S.element()][0]][0];
613
double b = begin_node()[begin()[S.element()][1]][0];
614
return point(a + S.hat()[0]*(b-a));
616
case geo_element::t: {
617
point a = begin_node()[begin()[S.element()][0]];
618
point b = begin_node()[begin()[S.element()][1]];
619
point c = begin_node()[begin()[S.element()][2]];
620
return ((!is_a_vector)?a:point(0,0)) + (b-a)*S.hat()[0]
623
case geo_element::T: {
624
point a = begin_node()[begin()[S.element()][0]];
625
point b = begin_node()[begin()[S.element()][1]];
626
point c = begin_node()[begin()[S.element()][2]];
627
point d = begin_node()[begin()[S.element()][3]];
628
return ((!is_a_vector)?a:point(0,0)) + (b-a)*S.hat()[0]
632
default : error_macro("dehatter: element not yet implemented: "
633
<< begin()[S.element()].name());
636
// ============================================================================
638
// ============================================================================
640
georep::localize (const point& x, geo_element::size_type& K_idx) const
642
switch (dimension()) {
644
K_idx = _h1d.localize(*this, x[0]);
645
return (K_idx != numeric_limits<hazel_1d::size_type>::max());
647
return _h.localize(x, K_idx);
649
warning_macro ("localize("<<x<<") = " << K_idx);
650
K_idx = _tree.search_element_containing (begin(), x);
651
return (K_idx != numeric_limits<size_t>::max());
652
default : error_macro("localizer not defined in " << dimension() << "D");
657
georep::localize_nearest (const point& x, point& y, geo_element::size_type& element) const
659
switch (dimension()) {
661
_h1d.localize_nearest (*this, x[0], y[0], element);
664
return _h.localize_nearest (x, y, element);
666
element = _tree.find_close_element (begin(), x, y);
668
default : error_macro("localizer/nearest not defined in " << dimension() << "D");
672
georep::trace (const point& x0, const point& v, point& x, Float& t, size_type& element) const
674
switch (dimension()) {
676
return _h1d.trace (*this, x0[0], v[0], x[0], t, element);
678
return _h.trace(x0, v, x, t, element);
680
default : error_macro("localizer/trace not defined in 3D");
684
// ============================================================================
685
// basic input/output
686
// ============================================================================
688
// ----------------------------------------------------------------------------
690
// ----------------------------------------------------------------------------
695
orheostream file (name(), "geo");
699
georep::dump (ostream& os) const
701
os << "dump mesh" << endl
702
<< "name " << _name << endl
703
<< "version " << _version << endl
704
<< "dimension " << _dim << endl
705
<< "xmin " << _xmin << endl
706
<< "xmax " << _xmax << endl
708
for (size_type i = 0; i < 4; i++)
709
os << "count_geo[" << i << "] = " << _count_geo[i] << endl;
710
for (size_type i = 0; i < reference_element::max_size; i++)
711
os << "count_element[" << i << "] = " << _count_element[i] << endl;
714
os << "n_vertex " << _x.size() << endl;
715
for (const_iterator_node iter_p = _x.begin(); iter_p != _x.end(); iter_p++, i++) {
716
os << "x[" << i << "] = ";
717
(*iter_p).put (os, _dim);
721
os << "n_element " << size() << endl;
722
for (const_iterator iter_e = begin(); iter_e != end(); iter_e++, i++) {
723
os << "E[" << i << "]: ";
726
os << "n_domain " << _domlist.size() << endl;
727
for (domlist_type::const_iterator i = _domlist.begin(); i != _domlist.end(); i++) {
732
ostream& georep::put(ostream& s) const
734
size_type map_d = map_dimension();
741
s << version() << " " << dimension() << " " << n_node() << " " << size();
742
if (version() >= 2) {
743
if (dimension() >= 3) {
745
s << " " << n_face();
750
if (dimension() >= 2) {
752
s << " " << n_edge();
759
if (version() >= 3) {
760
s << "coordinate_system " << coordinate_system() << endl;
766
if (_save_double_precision) s << setprecision(numeric_limits<Float>::digits10);
767
const_iterator_node iter_p = begin_node();
768
const_iterator_node last_p = end_node();
769
while (iter_p != last_p) {
770
const point& p = *iter_p++;
771
p.put (s, dimension());
778
const_iterator iter_e = begin();
779
const_iterator last_e = end();
780
while (iter_e != last_e) {
781
s << (*iter_e) << endl;
788
if (version() >= 2 && map_d >= 3) {
789
vector<geo_element> f (n_face());
790
const_iterator last_elt = end();
791
for (const_iterator iter_elt = begin(); iter_elt != last_elt; iter_elt++) {
792
const geo_element& K = *iter_elt;
793
for (size_type i = 0; i < K.n_face(); i++) {
794
size_type idx = K.subgeo (2, i);
795
if (f [idx].type() == reference_element::max_size) {
797
K.build_subgeo (2, i, f[idx]);
801
for (size_type j = 0; j < n_face(); j++) {
809
if (version() >= 2 && map_d >= 2) {
810
vector<size_type> v1(n_edge());
811
vector<size_type> v2(n_edge());
812
const_iterator last_elt = end();
813
for (const_iterator iter_elt = begin(); iter_elt != last_elt; iter_elt++) {
814
const geo_element& K = *iter_elt;
815
for (size_type i = 0; i < K.n_edge(); i++) {
816
v1 [K.edge(i)] = K.subgeo_vertex (1, i, 0);
817
v2 [K.edge(i)] = K.subgeo_vertex (1, i, 1);
820
for (size_type j = 0; j < n_edge(); j++) {
821
s << v1[j] << " " << v2[j] << "\n" ;
828
domlist_type::const_iterator iter = _domlist.begin();
829
domlist_type::const_iterator last = _domlist.end();
830
while (iter != last) s << endl << *iter++;
835
// =======================================================================
837
// =======================================================================
840
georep::check() const
842
warning_macro ("checking geo: elements...");
844
warning_macro ("checking geo done (version <= 1).");
847
for (const_iterator i = begin(); i != end(); i++) {
851
warning_macro ("checking geo: domains...");
852
domlist_type::const_iterator last = _domlist.end();
853
for (domlist_type::const_iterator i = _domlist.begin(); i != last; i++) {
855
for (domain::const_iterator j = (*i).begin(); j != (*i).end(); j++) {
860
warning_macro ("checking geo done.");
863
georep::check (const geo_element& E) const
865
for (size_type d = 0; d < E.dimension(); d++) {
866
for (size_type i_subgeo = 0; i_subgeo < E.n_subgeo(d); i_subgeo++) {
867
for (size_type i = 0; i < E.subgeo_size(d, i_subgeo); i++) {
868
size_type idx = E.subgeo_vertex(d, i_subgeo, i);
869
check_macro (idx < _count_geo[d], d << "D subgeo index `" << idx
870
<< " out of range 0:" << _count_geo[d]);
875
// ----------------------------------------------------------------------------
877
// ----------------------------------------------------------------------------
882
georep::iterator_node iter_p,
883
georep::iterator_node last_p,
886
georep::size_type dim)
888
typedef georep::size_type size_type;
889
for (size_type j = 0; j < dim; j++) {
890
xmin[j] = numeric_limits<Float>::max();
891
xmax[j] = -numeric_limits<Float>::max();
893
for (size_type j = dim+1; j < 3; j++) {
894
xmin [j] = xmax [j] = 0;
896
while (s.good() && iter_p != last_p) {
897
point& p = *iter_p++;
898
for (size_type j = 0 ; j < dim; j++) {
900
xmin[j] = min(p[j], xmin[j]);
901
xmax[j] = max(p[j], xmax[j]);
903
for (size_type j = dim; j < 3; j++)
909
georep::get_rheo(istream& is)
911
typedef georep::size_type size_type;
912
bool verbose = iorheo::getverbose(is);
914
if (!is.good()) error_macro("bad input stream for geo.");
916
if (!scatch(is,"\nmesh"))
917
error_macro("input stream does not contains a geo.");
926
is >> _version >> _dim >> n_vert >> n_elt;
930
warning_macro ("obsolete format version 1");
931
} else if (_version >= 2) {
939
fatal_macro("unknown mesh version " << _version);
942
if (!scatch(is,"\ncoordinate_system")) {
943
error_macro("read geo version 3: `coordinate_system' expected");
947
set_coordinate_system (sys_coord);
953
get_nodes (is, begin_node(), end_node(), _xmin, _xmax, _dim);
954
if (!is.good()) return is;
960
_count_geo [0] = n_node();
961
_count_element [0] = n_node();
962
georep::iterator last_elt = end();
963
for (georep::iterator i = begin(); is.good() && i != last_elt; i++, idx++) {
966
(*i).set_index (idx);
967
_count_geo [(*i).dimension()]++;
968
_count_element [(*i).type()]++;
970
size_type map_d = map_dimension();
972
vector<set<size_type> > ball (n_node());
973
if (_version >= 2 && map_d > 0) {
975
// TODO: _count_element [type]++
977
ball.resize(n_node());
978
build_point_to_element_sets (begin(), end(), ball.begin());
979
if (map_d >= 3) { // get faces
980
load_subgeo_numbering (begin(), ball.begin(), 2, n_fac, is,
981
_count_geo, _count_element);
983
load_subgeo_numbering (begin(), ball.begin(), 1, n_edg, is,
984
_count_geo, _count_element);
994
if (_version >= 2 && d.dimension() > 0 && size() > 0) {
995
propagate_subgeo_numbering (d.begin(), d.end(), begin(), ball.begin());
997
_domlist.push_back(d);
1002
operator >> (istream& s, georep& m)
1004
string basename = iorheo::getbasename(s);
1005
if (basename.length() != 0) {
1009
warning_macro ("unamed input mesh has been named `" << m._name << "'");
1011
m.reset_counters ();
1012
iorheo::flag_type fmt = iorheo::flags(s) & iorheo::format_field;
1013
if (fmt [iorheo::bamg]) { return m.get_bamg (s); }
1014
else if (fmt [iorheo::grummp]) { return m.get_grummp (s); }
1015
else if (fmt [iorheo::gmsh]) { return m.get_gmsh (s); }
1016
else if (fmt [iorheo::mmg3d]) { return m.get_mmg3d (s); }
1017
else if (fmt [iorheo::tetgen]) { return m.get_tetgen (s); }
1018
else if (fmt [iorheo::qmg]) { return m.get_qmg (s); }
1019
else if (fmt [iorheo::cemagref]) { return m.get_cemagref (s); }
1020
else if (fmt [iorheo::vtkpolydata]) { return m.get_vtk_polydata (s); }
1021
else { return m.get_rheo (s); }
1024
georep::may_have_version_2() const
1026
if (version() < 2) {
1027
error_macro("geo version 2 requiered (see: geo -upgrade): "
1028
<< _version << " founded");
1031
// ----------------------------------------------------------------------------
1033
// ----------------------------------------------------------------------------
1035
georep::erase_domain (const string& name)
1037
for (domlist_type::iterator q = _domlist.begin(); q != _domlist.end(); q++) {
1038
if (name == (*q).name()) {
1045
georep::insert_domain (const domain& d)
1047
_domlist.push_back(d);
1049
// ============================================================
1050
// construction of jump-interface domain data:
1051
// ============================================================
1053
georep::jump_interface(const domain& interface, const domain& region,
1054
map<size_type,tiny_vector<size_type> >& special_elements,
1055
map<size_type,size_type>& node_global_to_interface) const
1057
check_macro((_dim==2)&&(interface.dimension()==1), "Only 2D implemented");
1058
domain::const_iterator edge_i =interface.begin();
1059
domain::const_iterator end_edges =interface.end();
1060
set<size_type> n_tmp;
1061
for ( ; edge_i!=end_edges ; edge_i++)
1063
n_tmp.insert((*edge_i)[0]);
1064
n_tmp.insert((*edge_i)[1]);
1065
//node_pair* pair0 = new node_pair((*edge_i)[0], 0);
1066
//if (!(node_global_to_interface.insert(*pair0).second)) (*pair0).~pair();
1067
//node_pair* pair1 = new node_pair((*edge_i)[1], 0);
1068
//if (!(node_global_to_interface.insert(*pair1).second)) (*pair1).~pair();
1070
set<size_type>::iterator n_tmp_i =n_tmp.begin();
1071
set<size_type>::iterator end_n_tmp =n_tmp.end();
1072
for ( size_type j=0 ; n_tmp_i != end_n_tmp ; n_tmp_i++, j++ )
1074
//node_global_to_interface.insert(node_global_to_interface.end(),node_pair((*n_tmp_i),j));
1075
node_global_to_interface.insert(make_pair((*n_tmp_i),j));
1077
map<size_type,size_type>::iterator end_nodes =node_global_to_interface.end();
1078
map<size_type,tiny_vector<size_type> >::iterator special_K;
1079
domain::const_iterator ele_i =region.begin();
1080
domain::const_iterator end_elements =region.end();
1081
tiny_vector<size_type> *l_conn;
1082
for ( ; ele_i != end_elements ; ele_i++)
1084
bool already_inserted=false;
1085
for (size_type j=0; j <(*ele_i).size(); j++)
1087
map<size_type,size_type>::const_iterator node
1088
=node_global_to_interface.find((*ele_i)[j]);
1089
if (node!=end_nodes)
1091
if (!already_inserted)
1093
already_inserted=true;
1094
l_conn= new tiny_vector<size_type>((*ele_i).size());
1095
(*l_conn).fill(numeric_limits<size_type>::max());
1096
special_K = special_elements.insert(
1097
make_pair((*ele_i).index(),*l_conn)).first;
1099
(*special_K).second[j] = (*node).second;
1104
// ============================================================================
1106
// ============================================================================
1108
// small cpu-time equality check
1109
// we base on name changing convention
1110
// when mesh change !
1112
georep::operator == (const georep& y) const
1114
// rought connectivity comparison
1115
if (size() != y.size()) return false;
1116
if (n_node() != y.n_node()) return false;
1117
if (n_edge() != y.n_edge()) return false;
1118
if (n_face() != y.n_face()) return false;
1119
if (n_domain() != y.n_domain()) return false;
1121
// when only node coords change, we change the name !
1122
if (name() != y.name()) return false;
1127
// ============================================================================
1128
// Boundary builders
1129
// ============================================================================
1131
georep::boundary () const
1133
if (_boundary_initialized) return _boundary;
1134
_boundary_initialized = true;
1135
size_type d = dimension();
1136
_boundary.set_name("boundary");
1137
vector<set<size_type> > ball (n_node());
1138
build_point_to_element_sets (begin(), end(), ball.begin());
1139
build_boundary (begin(), end(), ball.begin(), d-1, _boundary);
1143
georep::label_interface(const domain& dom1, const domain& dom2, const string& name)
1145
// for 2D domains only
1146
check_macro(dom1.dimension()==2 && dom2.dimension()==2, "only 2D domains supported");
1150
interface.set_dimension(1);
1151
interface.set_name(name);
1152
// first step: sort vertices of larger domain in a set
1153
domain smaller, larger;
1154
if (dom1.size()<dom2.size()) { smaller=dom1; larger=dom2; }
1155
else { smaller=dom2; larger=dom1; }
1156
set< size_type > larger_vertices;
1157
domain::const_iterator iK=larger.begin();
1158
domain::const_iterator eK=larger.end();
1159
for ( ; iK!=eK; iK++)
1160
for (size_type i=0; i<(*iK).size(); i++)
1161
larger_vertices.insert((*iK)[i]);
1162
// second step: find smaller domain vertices which are within set and build interface
1165
for ( ; iK!=eK; iK++)
1166
for (size_type i=0, j=(*iK).size(); i<(*iK).size(); i++, j=i-1)
1167
if ( (larger_vertices.count((*iK)[j])>0)
1168
&& (larger_vertices.count((*iK)[i])>0) )
1170
K[0]=(*iK)[j]; K[1]=(*iK)[i];
1171
interface.push_back(K);
1173
_domlist.push_back(interface);
1176
georep::measure (const geo_element& K) const
1178
georep::const_iterator_vertex p = begin_vertex();
1180
case geo_element::p:
1182
case geo_element::e:
1183
switch (dimension()) {
1184
case 1: return fabs(p[K[1]][0] - p[K[0]][0]);
1185
case 2: return sqrt ( sqr(p[K[1]][0] - p[K[0]][0])
1186
+ sqr(p[K[1]][1] - p[K[0]][1]));
1187
default: return dist(p[K[1]], p[K[0]]);
1189
case geo_element::t:
1190
switch (dimension()) {
1191
case 2: return fabs(vect2d( p[K[1]] -p[K[0]], p[K[2]] -p[K[0]] )/2.0);
1192
default: return norm(vect ( p[K[1]] -p[K[0]], p[K[2]] -p[K[0]] )/2.0);
1194
case geo_element::q: {
1195
switch (dimension()) {
1197
Float area1 = fabs(vect2d( p[K[1]] -p[K[0]], p[K[2]] -p[K[0]] )/2);
1198
Float area2 = fabs(vect2d( p[K[2]] -p[K[0]], p[K[3]] -p[K[0]] )/2);
1199
return area1 + area2;
1202
Float area1 = norm(vect( p[K[1]] -p[K[0]], p[K[2]] -p[K[0]] )/2.0);
1203
Float area2 = norm(vect( p[K[2]] -p[K[0]], p[K[3]] -p[K[0]] )/2.0);
1204
return area1 + area2;
1208
case geo_element::T:
1209
return fabs(mixt(p[K[1]]-p[K[0]], p[K[2]]-p[K[0]], p[K[3]]-p[K[0]])/6);
1211
error_macro("area: element " << K.name() << " not yet supported yet.");
1215
// compute normal to a side of an element
1217
georep::normal (const geo_element& K, georep::size_type side) const
1220
case reference_element::e: { // works also in 2d and 3d
1221
point t = vertex(K[1]) - vertex(K[0]);
1223
if (side == 1) return t; else return -t;
1225
case reference_element::t:
1226
case reference_element::q: { // works also in 3d
1227
point t = vertex(K[(side+1)%3]) - vertex(K[side]);
1228
return point(t[1],-t[0])/norm(t);
1232
error_macro("normal to element not yet implemented for element type `"
1233
<< K.type() << "'");
1236
// compute normal to a facet: triangle in 3D, etc
1238
georep::normal (const geo_element& S) const
1241
case reference_element::e: { // 2d space
1242
point t = vertex(S[1]) - vertex(S[0]);
1244
return point (-t[1], t[0]);
1246
case reference_element::t:
1247
case reference_element::q: { // 3d space
1248
const point& a = vertex(S[0]);
1249
const point& b = vertex(S[1]);
1250
const point& c = vertex(S[2]);
1251
point n = vect (b-a, c-a);
1255
error_macro("normal to facet not yet implemented for element type `"
1256
<< S.type() << "'");
added
removed
nfem/lib/geo.cc
