1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
///
/// =========================================================================
// ref. ArnQin-1992 : for P2-P1d incompressible element
// split each triangle into 3 by inserting its barycenter
// example:
// mkgeo_grid -t 10 | ./geo_split_filter | geo -upgrade -geo - > isquare-10.geo
#include "rheolef.h"
using namespace rheolef;
using namespace std;
void split (const geo& omega) {
check_macro (omega.order() == 1, "only order 1 mesh still supported");
check_macro (omega.map_dimension() == 2, "invalid map_dimension");
size_t nn = omega.sizes().ownership_by_variant[reference_element::p].size();
size_t nt = omega.sizes().ownership_by_variant[reference_element::t].size();
check_macro (nt == omega.size(), "mesh may contain only triangles");
const disarray<point>& x = omega.get_nodes();
disarray<point> xg (nt);
for (size_t ie = 0; ie < nt; ++ie) {
const geo_element& K = omega[ie];
xg[ie] = point(0,0);
for (size_t iloc = 0; iloc < K.size(); ++iloc) {
xg[ie] += x[K[iloc]];
}
xg[ie] /= Float(int(K.size()));
}
cout << "#!geo" << endl
<< "mesh" << endl
<< "4" << endl
<< "header" << endl
<< " dimension 2" << endl
<< " nodes " << x.size() + nt << endl
<< " triangles " << 3*nt << endl;
if (omega.coordinate_system_name() != "cartesian") {
cout << " coordinate_system " << omega.coordinate_system_name() << endl;
}
cout << "end header" << endl
<< endl;
for (size_t i = 0; i < x.size(); ++i) {
x[i].put (cout,2);
cout << endl;
}
for (size_t i = 0; i < xg.size(); ++i) {
xg[i].put (cout,2);
cout << endl;
}
dout << endl;
for (size_t ie = 0; ie < nt; ++ie) {
const geo_element& K = omega[ie];
size_t ig = nn + ie;
for (size_t iloc = 0; iloc < K.size(); ++iloc) {
size_t iloc2 = (iloc+1) % K.size();
dout << "t\t" << K[iloc] << " " << K[iloc2] << " " << ig << endl;
}
}
// put domains: bdr edges are unchanged
for (size_t idom = 0; idom < omega.n_domain_indirect(); ++idom) {
const domain_indirect_basic<rheo_default_memory_model>& dom = omega.get_domain_indirect(idom);
dout << endl
<< "domain" << endl
<< dom.name() << endl
<< "1 1 " << dom.size() << endl;
for (size_t oige = 0; oige < dom.size(); ++oige) {
size_t ie = dom.oige(oige).index();
const geo_element& E = omega.get_geo_element(1,ie);
dout << "e\t" << E[0] << " " << E[1] << endl;
}
}
}
int main(int argc, char**argv) {
environment rheolef (argc, argv);
check_macro (communicator().size() == 1, "program may run sequentially");
geo omega;
din >> omega;
split(omega);
}
|
