geo_split.cc
Go to the documentation of this file.
// 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);
}
mkgeo_couette.nn
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:391
geo
see the geo page for the full documentation
split
void split(const geo &omega)
Definition: geo_split.cc:28