Rheolef  7.1
an efficient C++ finite element environment
geo_split.cc
Go to the documentation of this file.
1 // 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"
26 using namespace rheolef;
27 using namespace std;
28 void 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 }
87 int 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 }
mkgeo_couette.nn
int nn
Definition: mkgeo_couette.sh:71
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef.h
rheolef - reference manual
main
int main(int argc, char **argv)
Definition: geo_split.cc:87
rheolef::din
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::point
point_basic< Float > point
Definition: point.h:164
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::Float
double Float
see the Float page for the full documentation
Definition: Float.h:143
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
rheolef::domain_indirect_basic
the finite element boundary domain
Definition: domain_indirect.h:335
rheolef::std
Definition: vec_expr_v2.h:402
geo
see the geo page for the full documentation
split
void split(const geo &omega)
Definition: geo_split.cc:28