Loading...
Searching...
No Matches
geo_seq_get_bamg.cc
Go to the documentation of this file.
//
// geo: bamg input
//
// author: Pierre.Saramito@imag.fr
//
// 5 march 2012
//
#include "rheolef/geo.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
using namespace std;
namespace rheolef {
static
void
bamg_load_element (std::istream& is, size_t variant, geo_element_auto<>& K)
{
K.reset (variant, 1);
for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
is >> K[iloc];
K[iloc]--;
}
}
template <class T>
static
void
build_domains (
size_t dom_dim,
const disarray<geo_element_auto<>, sequential>& elt,
const disarray<size_t, sequential>& elt_bamg_dom_id,
const std::vector<std::string>& dom_name,
const geo_basic<T,sequential>& omega,
vector<index_set>* ball)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef pair<size_type,size_type> pair_t;
typedef geo_element_auto<> element_t;
typedef disarray<element_t, sequential> e_array_t;
typedef disarray<size_type, sequential> i_array_t;
typedef map<size_type, pair_t, less<size_type> > map_t;
// -------------------------------------------------------------------
// 1) reduce the bamg domain id to [0:ndom[ by counting domain id
// -------------------------------------------------------------------
// map[bamg_id] will gives idom and the number of its elements
map_t bamg_id2idom; // TODO (less<size_type>());
size_type idom = 0;
for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
size_type bamg_id = elt_bamg_dom_id [ie];
typename map_t::iterator iter = bamg_id2idom.find (bamg_id);
if (iter != bamg_id2idom.end()) {
// here is a previous bamg_dom_id: increment counter
((*iter).second.second)++;
continue;
}
// here is a new bamg_dom_id: associated to idom and elt counter=1
bamg_id2idom.insert (pair<size_type,pair_t>(bamg_id, pair_t(idom,1)));
idom++;
}
size_type ndom = bamg_id2idom.size();
// -------------------------------------------------------------------
// 2) check that ndom matches the domain name disarray size
// -------------------------------------------------------------------
if (ndom != dom_name.size()) {
warning_macro ("geo::get_bamg: "
<< dom_name.size() << " domain name(s) founded while "
<< ndom << " bamg "<<dom_dim<<"d domain(s) are defined");
cerr << endl;
error_macro("HINT: check domain name file (.dmn)");
}
// -------------------------------------------------------------------
// 3) create domain disarray: loop on elements
// -------------------------------------------------------------------
element_t dummy_elt (reference_element::p, 0);
vector<e_array_t> domain (ndom, e_array_t(0, dummy_elt));
vector<size_type> counter (ndom, 0);
for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
size_type bamg_dom_id = elt_bamg_dom_id [ie];
size_type idom = bamg_id2idom [bamg_dom_id].first;
if (domain[idom].size() == 0) {
size_type size = bamg_id2idom [bamg_dom_id].second;
domain[idom] = e_array_t (size, dummy_elt);
}
size_type dom_ie = counter[idom];
domain[idom][dom_ie] = elt[ie];
counter[idom]++;
}
// -------------------------------------------------------------------
// 3) insert domains in the mesh: loop on domains
// -------------------------------------------------------------------
for (typename map_t::const_iterator
iter = bamg_id2idom.begin(),
last = bamg_id2idom.end(); iter != last; ++iter) {
size_type bamg_id = (*iter).first;
size_type idom = (*iter).second.first;
size_type size = (*iter).second.second;
check_macro (idom < dom_name.size(), "invalid idom="<<idom<<" for domain name");
string name = dom_name [idom];
domain_indirect_basic<sequential> dom (domain[idom], omega, ball);
dom.set_name (name);
dom.set_map_dimension (dom_dim);
omega.insert_domain_indirect (dom);
}
size_type ndom2 = omega.n_domain_indirect();
}
template <class T>
static
void
build_vertex_domains (
const disarray<size_t, sequential>& edg_bdr_bamg_dom_id,
const disarray<size_t, sequential>& vert_bamg_dom_id,
const std::vector<std::string>& dom_name,
const geo_basic<T,sequential>& omega,
vector<index_set>* ball)
{
if (dom_name.size() == 0) return;
typedef geo_element_auto<> element_t;
typedef disarray<element_t, sequential> v_array_t;
element_t dummy_elt (reference_element::p, 0);
// -------------------------------------------------------------------
// 1) list all vertex domain id
// -------------------------------------------------------------------
std::set<size_t> vert_id;
for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
size_t dom_id = vert_bamg_dom_id [iv];
if (dom_id == 0) continue;
vert_id.insert (dom_id);
}
// -------------------------------------------------------------------
// 2) omit vertex that are marked with an edge id
// -------------------------------------------------------------------
for (size_t iedg_bdr = 0, nedg_bdr = edg_bdr_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
size_t dom_id = edg_bdr_bamg_dom_id [iedg_bdr];
vert_id.erase (dom_id);
}
check_macro (vert_id.size() == dom_name.size(),
"unexpected VertexDomainNames with "<<dom_name.size()
<<" domain names while the mesh provides " << vert_id.size()
<<" vertex labels");
// -------------------------------------------------------------------
// 3) loop on vertex domain and insert it in mesh
// -------------------------------------------------------------------
size_t idom = 0;
for (std::set<size_t>::const_iterator iter = vert_id.begin(); iter != vert_id.end(); ++iter, ++idom) {
size_t id = *iter;
string name = dom_name[idom];
size_t dom_size = 0;
for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
if (vert_bamg_dom_id [iv] == id) dom_size++;
}
v_array_t vert_list (dom_size, dummy_elt);
for (size_t iv = 0, iv_dom = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
if (vert_bamg_dom_id [iv] != id) continue;
element_t& K = vert_list [iv_dom];
K.reset (reference_element::p, 1);
K[0] = iv;
iv_dom++;
}
domain_indirect_basic<sequential> dom (vert_list, omega, ball);
dom.set_name (name);
dom.set_map_dimension (0);
omega.insert_domain_indirect (dom);
}
}
template <class T>
idiststream&
geo_get_bamg (idiststream& ips, geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef typename geo_basic<T,sequential>::node_type node_type;
check_macro (ips.good(), "bad input stream for bamg");
istream& is = ips.is();
geo_header hdr;
hdr.dimension = 2;
hdr.map_dimension = 2;
hdr.order = 1;
hdr.dis_size_by_dimension [2] = 0;
// ------------------------------------------------------------------------
// 1) load data
// ------------------------------------------------------------------------
typedef geo_element_auto<> element_t;
typedef disarray<node_type, sequential> n_array_t;
typedef disarray<element_t, sequential> e_array_t;
typedef disarray<size_type, sequential> i_array_t;
n_array_t vertex;
i_array_t vert_bamg_dom_id;
i_array_t edg_bdr_bamg_dom_id;
i_array_t tri_bamg_dom_id;
i_array_t qua_bamg_dom_id;
e_array_t edg_bdr;
std::array<e_array_t, reference_element::max_variant> tmp_geo_element;
string label;
while (is.good() && label != "End") {
is >> label;
if (label == "Vertices") {
// ------------------------------------------------------------------------
// 1.1) load the coordinates
// Vertices <np>
// {xi yi dom_idx} i=0..np-1
// ------------------------------------------------------------------------
size_type nnod;
is >> nnod;
if (nnod == 0) {
warning_macro("empty bamg mesh file");
return ips;
}
hdr.dis_size_by_dimension [0] = nnod;
hdr.dis_size_by_variant [0] = nnod;
vertex.resize (nnod);
vert_bamg_dom_id = i_array_t (nnod, 0);
for (size_type inod = 0; inod < nnod; inod++) {
is >> vertex[inod][0] >> vertex[inod][1] >> vert_bamg_dom_id[inod];
}
check_macro (ips.good(), "bad input stream for bamg");
} else if (label == "Triangles") {
// ------------------------------------------------------------------------
// 2.1) load the connectivity
// Triangle <nt>
// {s0 s1 s2 dom_idx} j=0..nt-1
// ------------------------------------------------------------------------
size_type nt;
is >> nt;
hdr.dis_size_by_dimension [2] += nt;
hdr.dis_size_by_variant [reference_element::t] = nt;
element_t init_tri (reference_element::t, hdr.order);
tmp_geo_element [reference_element::t] = e_array_t (nt, init_tri);
tri_bamg_dom_id = i_array_t (nt, 0);
for (size_type it = 0; it < nt; it++) {
bamg_load_element (is, reference_element::t,
tmp_geo_element[reference_element::t][it]);
is >> tri_bamg_dom_id[it];
}
} else if (label == "Quadrilaterals") {
// ------------------------------------------------------------------------
// Quadrilaterals <nq>
// {s0 s1 s2 s3 dom_idx} j=0..nq-1
// ------------------------------------------------------------------------
size_type nq;
is >> nq;
hdr.dis_size_by_dimension [2] += nq;
hdr.dis_size_by_variant [reference_element::q] = nq;
element_t init_qua (reference_element::q, hdr.order);
tmp_geo_element [reference_element::q] = e_array_t (nq, init_qua);
qua_bamg_dom_id = i_array_t (nq, 0);
for (size_type iq = 0; iq < nq; iq++) {
bamg_load_element (is, reference_element::q,
tmp_geo_element[reference_element::q][iq]);
is >> qua_bamg_dom_id[iq];
}
} else if (label == "Edges") {
// ------------------------------------------------------------------------
// 2.3) load the boundary domains
// Edges <nedg>
// {s0 s1 dom_idx} j=0..nedg-1
// ------------------------------------------------------------------------
size_type nedg;
is >> nedg;
element_t init_edg (reference_element::e, hdr.order);
edg_bdr = e_array_t (nedg, init_edg);
edg_bdr_bamg_dom_id = i_array_t (nedg, 0);
for (size_type iedg = 0; iedg < nedg; iedg++) {
element_t& K = edg_bdr[iedg];
K.reset (reference_element::e, hdr.order);
for (size_type iloc = 0, nloc = 2; iloc < nloc; iloc++) {
is >> K[iloc];
K[iloc]--;
}
is >> edg_bdr_bamg_dom_id[iedg];
}
}
}
// ---------------------------------------------------------------
// 2) check rheolef extension to optional domain names
// ---------------------------------------------------------------
vector<string> vertice_domain_name;
vector<string> edg_dom_name;
vector<string> region_domain_name;
char c;
is >> ws >> c; // skip white and grab a char
// have "EdgeDomainNames" or "VertexDomainNames" ?
// bamg mesh may be followed by field data and such, so be carrefull...
while (c == 'E' || c == 'V' || c == 'R') {
is.unget(); // put char back
if (c == 'R') {
if (!scatch(is,"RegionDomainNames")) break;
size_type n_dom_region;
is >> n_dom_region;
region_domain_name.resize (n_dom_region);
for (size_type k = 0; k < n_dom_region; k++) {
is >> region_domain_name[k];
}
} else if (c == 'E') {
if (!scatch(is,"EdgeDomainNames")) break;
// ---------------------------------------------------------------
// get edge domain names
// ---------------------------------------------------------------
size_type n_dom_edge;
is >> n_dom_edge;
edg_dom_name.resize (n_dom_edge);
for (size_type k = 0; k < n_dom_edge; k++) {
is >> edg_dom_name[k];
}
} else if (c == 'V') {
if (!scatch(is,"VertexDomainNames")) break;
// ---------------------------------------------------------------
// get vertice domain names
// ---------------------------------------------------------------
size_type n_dom_vertice;
is >> n_dom_vertice;
vertice_domain_name.resize (n_dom_vertice);
for (size_type k = 0; k < n_dom_vertice; k++) {
is >> vertice_domain_name[k];
}
}
is >> ws >> c; // skip white and grab a char
}
// ------------------------------------------------------------------------
// 3) build & upgrade
// ------------------------------------------------------------------------
bool do_upgrade = true;
omega.build_from_data (hdr, vertex, tmp_geo_element, do_upgrade);
// ------------------------------------------------------------------------
// 4) get domain, until end-of-file
// ------------------------------------------------------------------------
vector<index_set> ball[4];
build_domains (1, edg_bdr, edg_bdr_bamg_dom_id, edg_dom_name, omega, ball);
build_vertex_domains (edg_bdr_bamg_dom_id, vert_bamg_dom_id, vertice_domain_name, omega, ball);
//TODO: region(d=2) domains
return ips;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template idiststream& geo_get_bamg<Float> (idiststream&, geo_basic<Float,sequential>&);
}// namespace rheolef
size_type
field::size_type size_type
Definition branch.cc:430
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element_auto
Definition geo_element.h:376
rheolef::idiststream
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
rheolef::idiststream::is
std::istream & is()
Definition diststream.h:400
rheolef::idiststream::good
bool good() const
Definition diststream.cc:124
rheolef::reference_element::q
static const variant_type q
Definition reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition reference_element.h:76
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::domain
geo domain
Definition geo_domain_indirect.h:58
error_macro
#define error_macro(message)
Definition dis_macros.h:49
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
vertex
const point vertex[n_vertex]
Definition edge.icc:67
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_contraction.name
str name
Definition mkgeo_contraction.sh:133
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::geo_get_bamg
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
Definition geo_seq_get_bamg.cc:189
rheolef::scatch
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition scatch.icc:44
rheolef::geo_get_bamg< Float >
template idiststream & geo_get_bamg< Float >(idiststream &, geo_basic< Float, sequential > &)
std
STL namespace.
rheolef::geo_header
Definition geo_header.h:32
rheolef::geo_header::map_dimension
size_type map_dimension
Definition geo_header.h:40
rheolef::geo_header::dis_size_by_dimension
size_type dis_size_by_dimension[4]
Definition geo_header.h:44
rheolef::geo_header::order
size_type order
Definition geo_header.h:42
rheolef::geo_header::dimension
size_type dimension
Definition geo_header.h:39
rheolef::geo_header::dis_size_by_variant
size_type dis_size_by_variant[reference_element::max_variant]
Definition geo_header.h:43