Loading...
Searching...
No Matches
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);
}
Float
see the Float page for the full documentation
communicator
see the communicator page for the full documentation
geo
see the geo page for the full documentation
point
see the point page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::domain_indirect_basic
the finite element boundary domain
Definition domain_indirect.h:335
rheolef::environment
see the environment page for the full documentation
Definition environment.h:121
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_element::size
size_type size() const
Definition geo_element.h:168
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
main
int main()
Definition field2bb.cc:58
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
split
void split(const geo &omega)
Definition geo_split.cc:28
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
std
STL namespace.
rheolef.h
rheolef - reference manual