Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_split.cc
Go to the documentation of this file.
1
21// ref. ArnQin-1992 : for P2-P1d incompressible element
22// split each triangle into 3 by inserting its barycenter
23// example:
24// mkgeo_grid -t 10 | ./geo_split_filter | geo -upgrade -geo - > isquare-10.geo
25#include "rheolef.h"
26using namespace rheolef;
27using namespace std;
28void split (const geo& omega) {
29 check_macro (omega.order() == 1, "only order 1 mesh still supported");
30 check_macro (omega.map_dimension() == 2, "invalid map_dimension");
31 size_t nn = omega.sizes().ownership_by_variant[reference_element::p].size();
32 size_t nt = omega.sizes().ownership_by_variant[reference_element::t].size();
33 check_macro (nt == omega.size(), "mesh may contain only triangles");
34 const disarray<point>& x = omega.get_nodes();
35 disarray<point> xg (nt);
36 for (size_t ie = 0; ie < nt; ++ie) {
37 const geo_element& K = omega[ie];
38 xg[ie] = point(0,0);
39 for (size_t iloc = 0; iloc < K.size(); ++iloc) {
40 xg[ie] += x[K[iloc]];
41 }
42 xg[ie] /= Float(int(K.size()));
43 }
44 cout << "#!geo" << endl
45 << "mesh" << endl
46 << "4" << endl
47 << "header" << endl
48 << " dimension 2" << endl
49 << " nodes " << x.size() + nt << endl
50 << " triangles " << 3*nt << endl;
51 if (omega.coordinate_system_name() != "cartesian") {
52 cout << " coordinate_system " << omega.coordinate_system_name() << endl;
53 }
54 cout << "end header" << endl
55 << endl;
56 for (size_t i = 0; i < x.size(); ++i) {
57 x[i].put (cout,2);
58 cout << endl;
59 }
60 for (size_t i = 0; i < xg.size(); ++i) {
61 xg[i].put (cout,2);
62 cout << endl;
63 }
64 dout << endl;
65 for (size_t ie = 0; ie < nt; ++ie) {
66 const geo_element& K = omega[ie];
67 size_t ig = nn + ie;
68 for (size_t iloc = 0; iloc < K.size(); ++iloc) {
69 size_t iloc2 = (iloc+1) % K.size();
70 dout << "t\t" << K[iloc] << " " << K[iloc2] << " " << ig << endl;
71 }
72 }
73 // put domains: bdr edges are unchanged
74 for (size_t idom = 0; idom < omega.n_domain_indirect(); ++idom) {
75 const domain_indirect_basic<rheo_default_memory_model>& dom = omega.get_domain_indirect(idom);
76 dout << endl
77 << "domain" << endl
78 << dom.name() << endl
79 << "1 1 " << dom.size() << endl;
80 for (size_t oige = 0; oige < dom.size(); ++oige) {
81 size_t ie = dom.oige(oige).index();
82 const geo_element& E = omega.get_geo_element(1,ie);
83 dout << "e\t" << E[0] << " " << E[1] << endl;
84 }
85 }
86}
87int main(int argc, char**argv) {
88 environment rheolef (argc, argv);
89 check_macro (communicator().size() == 1, "program may run sequentially");
90 geo omega;
91 din >> omega;
92 split(omega);
93}
see the Float page for the full documentation
see the communicator page for the full documentation
see the geo page for the full documentation
see the point page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
the finite element boundary domain
see the environment page for the full documentation
see the geo_element page for the full documentation
size_type size() const
static const variant_type p
static const variant_type t
int main()
Definition field2bb.cc:58
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void split(const geo &omega)
Definition geo_split.cc:28
This file is part of Rheolef.
STL namespace.
rheolef - reference manual