Loading...
Searching...
No Matches
geo_seq_get.cc
Go to the documentation of this file.
#include "rheolef/geo.h"
#include "rheolef/geo_domain.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/rheostream.h"
namespace rheolef {
template <class T> idiststream& geo_get_vtk (idiststream& ips, geo_basic<T,sequential>& omega);
template <class T> idiststream& geo_get_bamg (idiststream& ips, geo_basic<T,sequential>& omega);
template <class T>
idiststream&
geo_basic<T,sequential>::get (idiststream& ips)
{
iorheo::flag_type format = iorheo::flags(ips.is()) & iorheo::format_field;
if (format [iorheo::vtk]) { return geo_get_vtk (ips,*this); }
if (format [iorheo::bamg]) { return geo_get_bamg (ips,*this); }
// else: standard .geo format:
// allocate a new geo_rep object (TODO: do a dynamic_cast ?)
geo_rep<T,sequential>* ptr = new_macro((geo_rep<T,sequential>));
ptr->get (ips);
base::operator= (ptr);
return ips;
}
template <class T>
void
geo_rep<T,sequential>::set_element_side_index (size_type side_dim)
{
if (map_dimension() <= side_dim) return;
// ------------------------------------------------------------------------
// 1) ball(X) := { E; X is a vertex of E }
// ------------------------------------------------------------------------
index_set empty_set; // TODO: add a global allocator to empty_set
disarray<index_set,sequential> ball (geo_element_ownership(0), empty_set);
{
size_type isid = 0;
for (const_iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
const geo_element& S = *iter;
for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
size_type iv = S[iloc];
ball [iv] += isid;
}
}
}
// ------------------------------------------------------------------------
// 2) pour K dans partition(iproc)
// pour (dis_A,dis_B) arete de K
// set = dis_ball(dis_A) inter dis_ball(dis_B) = {dis_iedg}
// E = dis_edges(dis_iedg)
// => on numerote dis_iedg cette arete dans le geo_element K
// et on indique son orient en comparant a E, arete qui definit l'orient
// ------------------------------------------------------------------------
for (size_type dim = side_dim+1; dim <= base::_gs._map_dimension; dim++) {
for (iterator iter = begin(dim), last = end(dim); iter != last; iter++) {
geo_element& K = *iter;
for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
size_type jv0 = K[loc_jv0];
index_set isid_set = ball [jv0]; // copy: will be intersected
for (size_type sid_jloc = 1, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
size_type jv = K[loc_jv];
const index_set& ball_jv = ball [jv];
isid_set.inplace_intersection (ball_jv);
}
check_macro (isid_set.size() == 1, "connectivity: side not found in the side set");
size_type isid = *(isid_set.begin());
const geo_element& S = get_geo_element(side_dim,isid);
if (side_dim == 1) {
// side: edge
size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
geo_element::orientation_type orient = S.get_edge_orientation (jv0, jv1);
K.edge_indirect (loc_isid).set (orient, isid);
} else { // side_dim == 2
geo_element::orientation_type orient;
geo_element::shift_type shift;
if (K.subgeo_size (side_dim, loc_isid) == 3) {
// side: triangle
size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
S.get_orientation_and_shift (jv0, jv1, jv2, orient, shift);
} else {
// side: quadrangle
size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
size_type jv3 = K [K.subgeo_local_vertex (side_dim, loc_isid, 3)];
S.get_orientation_and_shift (jv0, jv1, jv2, jv3, orient, shift);
}
K.face_indirect (loc_isid).set (orient, isid, shift);
}
}
}
}
}
// =========================================================================
// get
// =========================================================================
template <class T>
idiststream&
geo_rep<T,sequential>::get (idiststream& ips)
{
using namespace std;
check_macro (ips.good(), "bad input stream for geo.");
istream& is = ips.is();
// ------------------------------------------------------------------------
// 0) get header
// ------------------------------------------------------------------------
check_macro (dis_scatch(ips,ips.comm(),"\nmesh"), "input stream does not contains a geo.");
ips >> base::_version;
check_macro (base::_version == 4, "mesh format version 4 expected, but format version " << base::_version << " founded");
geo_header hdr;
ips >> hdr;
bool do_upgrade = iorheo::getupgrade(is);
if (do_upgrade || hdr.need_upgrade()) {
return get_upgrade (ips, hdr);
} else {
return get_standard (ips, hdr);
}
}
template <class T>
idiststream&
geo_rep<T,sequential>::get_standard (idiststream& ips, const geo_header& hdr)
{
using namespace std;
check_macro (ips.good(), "bad input stream for geo.");
istream& is = ips.is();
// ------------------------------------------------------------------------
// 1) store header infos in geo
// ------------------------------------------------------------------------
base::_have_connectivity = true;
base::_name = "unnamed";
base::_dimension = hdr.dimension;
base::_gs._map_dimension = hdr.map_dimension;
base::_sys_coord = hdr.sys_coord;
base::_piola_basis.reset_family_index (hdr.order);
size_type nnod = hdr.dis_size_by_dimension [0];
size_type n_edg = hdr.dis_size_by_dimension [1];
size_type n_fac = hdr.dis_size_by_dimension [2];
size_type n_elt = hdr.dis_size_by_dimension [base::_gs._map_dimension];
// ------------------------------------------------------------------------
// 2) get coordinates
// ------------------------------------------------------------------------
base::_node.resize (nnod);
if (base::_dimension > 0) {
base::_node.get_values (ips, _point_get<T>(geo_base_rep<T,sequential>::_dimension));
check_macro (ips.good(), "bad input stream for geo.");
}
base::_gs.node_ownership = base::_node.ownership();
// ------------------------------------------------------------------------
// 3) get elements
// ------------------------------------------------------------------------
if (base::_gs._map_dimension > 0) {
for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
geo_element::parameter_type param (variant, 1);
base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
base::_geo_element [variant].get_values (ips);
base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
}
base::_gs.ownership_by_dimension [base::_gs._map_dimension] = distributor (n_elt, base::comm(), n_elt);
}
// ------------------------------------------------------------------------
// 4) check that nodes are numbered by increasing node_subgeo_dim
// ------------------------------------------------------------------------
// ICI va devenir obsolete car les noeuds seront numerotes par _numbering=Pk_numbering
{
std::vector<size_type> node_subgeo_dim (nnod, size_type(-1));
size_type prev_variant = 0;
for (iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++) {
geo_element& K = *iter;
check_macro (prev_variant <= K.variant(), "elements should be numbered by increasing variants (petqTPH)");
prev_variant = K.variant();
for (size_type d = 0; d <= base::_gs._map_dimension; d++) {
for (size_type loc_inod = K.first_inod(d), loc_nnod = K.last_inod(d); loc_inod < loc_nnod; loc_inod++) {
node_subgeo_dim [K[loc_inod]] = d;
}
}
}
size_type prev_node_dim = 0;
for (typename std::vector<size_type>::const_iterator iter = node_subgeo_dim.begin(), last = node_subgeo_dim.end();
iter != last; iter++) {
check_macro (prev_node_dim <= *iter, "nodes should be numbered by increasing subgeo dimension");
prev_node_dim = *iter;
}
}
// ------------------------------------------------------------------------
// 5) compute n_vert (n_vert < nnod when order > 1) and set element indexes (K.dis_ie & K.ios_dis_ie)
// ------------------------------------------------------------------------
size_type n_vert = 0;
if (base::_gs._map_dimension == 0) {
n_vert = nnod;
} else {
std::vector<size_t> is_vertex (nnod, 0);
size_type ie = 0;
for (iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++, ie++) {
geo_element& K = *iter;
K.set_ios_dis_ie (ie);
K.set_dis_ie (ie);
if (base::order() > 1) {
for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
is_vertex [K[iloc]] = 1;
}
}
}
if (base::order() == 1) {
n_vert = nnod;
} else {
n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
}
}
// ------------------------------------------------------------------------
// 6) create vertex-element (0d elements)
// ------------------------------------------------------------------------
geo_element::parameter_type param (reference_element::p, 1);
base::_geo_element [reference_element::p].resize (n_vert, param);
size_type first_iv = base::_node.ownership().first_index();
{
size_type iv = 0;
for (iterator iter = begin(0), last = end(0); iter != last; iter++, iv++) {
geo_element& P = *iter;
P[0] = first_iv + iv;
P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
P.set_ios_dis_ie (first_iv + iv);
}
}
// ownership_by_dimension[0]: used by connectivity
base::_gs.ownership_by_variant [reference_element::p] = base::_geo_element [reference_element::p].ownership();
base::_gs.ownership_by_dimension [0] = base::_geo_element [reference_element::p].ownership();
// ------------------------------------------------------------------------
// 7) get faces & edge
// ------------------------------------------------------------------------
if (base::_gs._map_dimension > 0) {
for (size_type side_dim = base::_gs._map_dimension - 1; side_dim >= 1; side_dim--) {
for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
geo_element::parameter_type param (variant, 1);
base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
base::_geo_element [variant].get_values (ips);
base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
}
size_type nsid = hdr.dis_size_by_dimension [side_dim];
base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
size_type isid = 0;
for (iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
geo_element& S = *iter;
S.set_ios_dis_ie (isid);
S.set_dis_ie (isid);
}
}
}
// ------------------------------------------------------------------------
// 8) get domain, until end-of-file
// ------------------------------------------------------------------------
vector<index_set> ball [4];
domain_indirect_basic<sequential> dom;
while (dom.get (ips, *this, ball)) {
base::_domains.push_back (dom);
}
// ------------------------------------------------------------------------
// 9) set indexes on faces and edges of elements, for P2 approx
// ------------------------------------------------------------------------
set_element_side_index (1);
set_element_side_index (2);
// ------------------------------------------------------------------------
// 10) bounding box: _xmin, _xmax
// ------------------------------------------------------------------------
base::compute_bbox();
return ips;
}
// ----------------------------------------------------------------------------
// dump
// ----------------------------------------------------------------------------
template <class T>
void
geo_rep<T,sequential>::dump (std::string name) const {
std::ofstream os ((name + "-dump.geo").c_str());
odiststream ods (os, base::_node.ownership().comm());
put_geo (ods);
}
// ----------------------------------------------------------------------------
// read from file
// ----------------------------------------------------------------------------
template <class T>
void
geo_rep<T,sequential>::load (std::string filename, const communicator&)
{
idiststream ips;
ips.open (filename, "geo");
check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
get (ips);
std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
std::string name = get_basename (root_name);
base::_name = name;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_rep<Float,sequential>;
template idiststream& geo_basic<Float,sequential>::get (idiststream&);
} // namespace rheolef
size_type
field::size_type size_type
Definition branch.cc:430
communicator
see the communicator page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::domain_indirect_basic
the finite element boundary domain
Definition domain_indirect.h:335
rheolef::geo_base_rep
base class for M=sequential or distributed meshes representations
Definition geo.h:528
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element_indirect::set
void set(orientation_type orient, size_type ige, size_type shift=0)
Definition geo_element_indirect.h:69
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition geo_element.h:134
rheolef::geo_element::get_orientation_and_shift
bool get_orientation_and_shift(const geo_element &S, orientation_type &orient, shift_type &shift) const
return orientation and shift between *this element and S
Definition geo_element.cc:114
rheolef::geo_element::subgeo_local_vertex
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
Definition geo_element.h:223
rheolef::geo_element::first_inod
size_type first_inod(size_type subgeo_dim) const
Definition geo_element.h:225
rheolef::geo_element::size
size_type size() const
Definition geo_element.h:168
rheolef::geo_element::set_ios_dis_ie
void set_ios_dis_ie(size_type ios_dis_ie)
Definition geo_element.h:173
rheolef::geo_element::subgeo_size
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
Definition geo_element.h:221
rheolef::geo_element::edge_indirect
const geo_element_indirect & edge_indirect(size_type i) const
Definition geo_element.h:193
rheolef::geo_element::variant
variant_type variant() const
Definition geo_element.h:161
rheolef::geo_element::get_edge_orientation
orientation_type get_edge_orientation(size_type dis_iv0, size_type dis_iv1) const
Definition geo_element.cc:155
rheolef::geo_element::shift_type
geo_element_indirect::shift_type shift_type
Definition geo_element.h:135
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition geo_element.h:212
rheolef::geo_element::face_indirect
const geo_element_indirect & face_indirect(size_type i) const
Definition geo_element.h:197
rheolef::geo_element::last_inod
size_type last_inod(size_type subgeo_dim) const
Definition geo_element.h:227
rheolef::geo_element::set_dis_ie
void set_dis_ie(size_type dis_ie)
Definition geo_element.h:172
rheolef::geo_rep< T, sequential >::iterator
base::iterator iterator
Definition geo.h:796
rheolef::geo_rep< T, sequential >::size_type
base::size_type size_type
Definition geo.h:791
rheolef::geo_rep< T, sequential >::const_iterator
base::const_iterator const_iterator
Definition geo.h:797
rheolef::geo_rep
sequential mesh representation
Definition geo.h:778
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::open
void open(std::string filename, std::string suffix="", const communicator &comm=communicator())
This routine opens a physical input file.
Definition diststream.cc:85
rheolef::idiststream::good
bool good() const
Definition diststream.cc:124
rheolef::idiststream::comm
const communicator & comm() const
Definition diststream.h:356
rheolef::index_set
Definition index_set_header.icc:58
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition index_set_body.icc:126
rheolef::iorheo::flag_type
std::bitset< last > flag_type
Definition iorheo.h:440
rheolef::iorheo::format_field
static flag_type format_field
Definition iorheo.h:445
rheolef::iorheo::vtk
@ vtk
Definition iorheo.h:408
rheolef::iorheo::bamg
@ bamg
Definition iorheo.h:390
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition reference_element.h:150
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition reference_element.h:148
rheolef::reference_element::p
static const variant_type p
Definition reference_element.h:75
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::delete_suffix
string delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
Definition rheostream.cc:226
rheolef::geo_get_bamg
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
Definition geo_seq_get_bamg.cc:189
rheolef::get_basename
string get_basename(const string &name)
get_basename: see the rheostream page for the full documentation
Definition rheostream.cc:258
rheolef::geo_get_vtk
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
Definition geo_seq_get_vtk.cc:64
rheolef::dis_scatch
bool dis_scatch(idiststream &ips, const communicator &comm, std::string ch)
distributed version of scatch(istream&,string)
Definition diststream.cc:44
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
rheolef::_point_get
point input helper
Definition geo.h:155
rheolef::geo_element::parameter_type
Definition geo_element.h:136
rheolef::geo_header
Definition geo_header.h:32
rheolef::geo_header::need_upgrade
bool need_upgrade() const
Definition geo_header.cc:79
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::sys_coord
coordinate_type sys_coord
Definition geo_header.h:41
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