21#include "rheolef/geo.h"
22#include "rheolef/geo_domain.h"
23#include "rheolef/space_numbering.h"
24#include "rheolef/rheostream.h"
25#include "rheolef/iorheo.h"
35template <
class T,
class M>
44 _have_connectivity(false),
45 _have_neighbour(false),
48 _sys_coord(space_constant::cartesian),
55 _tracer_ray_boundary(),
59template <
class T,
class M>
63 _version (o._version),
64 _serial_number (o._serial_number),
65 _geo_element (o._geo_element),
67 _domains (o._domains),
68 _have_connectivity (o._have_connectivity),
69 _have_neighbour (o._have_neighbour),
71 _dimension (o._dimension),
72 _sys_coord (o._sys_coord),
77 _piola_basis (o._piola_basis),
78 _locator (o._locator),
79 _tracer_ray_boundary (o._tracer_ray_boundary),
80 _nearestor (o._nearestor)
82 trace_macro (
"*** PHYSICAL COPY OF GEO_BASE_REP ***");
93 trace_macro (
"*** PHYSICAL COPY OF GEO_REP<seq> ***");
101 return new_macro (rep(*
this));
110#ifdef _RHEOLEF_HAVE_MPI
114 _inod2ios_dis_inod(),
115 _ios_inod2dis_inod(),
118 _igev2ios_dis_igev(),
125 _inod2ios_dis_inod(o._inod2ios_dis_inod),
126 _ios_inod2dis_inod(o._ios_inod2dis_inod),
127 _ios_ige2dis_ige (o._ios_ige2dis_ige),
129 _igev2ios_dis_igev(o._igev2ios_dis_igev),
130 _ios_igev2dis_igev(o._ios_igev2dis_igev)
132 trace_macro (
"*** PHYSICAL COPY OF GEO_REP<dis> ***");
140 return new_macro (rep(*
this));
146template<
class T,
class M>
149template<
class T,
class M>
161is_domain (
const std::string& name, std::string& bgd_name, std::string& dom_name)
163 size_t n = name.length();
164 if (n <= 2 || name[n-1] !=
']')
return false;
166 for (; i > 0 && name[i] !=
'[' && name[i] !=
'/'; i--) ;
168 if (i == 0 || i == n-2 || name[i] ==
'/')
return false;
170 bgd_name = name.substr(0,i);
171 dom_name = name.substr(i+1);
172 dom_name = dom_name.substr(0,dom_name.length()-1);
175template<
class T,
class M>
181 std::string bgd_name, dom_name;
182 if (is_domain(root_name,bgd_name,dom_name)) {
184 return omega[dom_name];
192#ifdef DEBUG_LOADED_TABLE_TRACE
196 omega.base::operator= (base((*iter).second,
typename base::internal()));
203 omega.geo_basic<
T,
M>::base::operator= (ptr);
206#ifdef DEBUG_LOADED_TABLE_TRACE
217template <
class T,
class M>
229 new_gamma.geo_basic<
T,
M>::base::operator= (geo_ptr);
238#define _RHEOLEF_save(M) \
241geo_basic<T,M>::save (std::string filename) const \
243 if (filename == "") filename = name(); \
244 odiststream out (filename, "geo"); \
247#define _RHEOLEF_set_name(M) \
250geo_basic<T,M>::set_name (std::string new_name) \
252 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
253 check_macro (ptr != 0, "cannot set_name on geo_domains"); \
254 ptr->set_name(new_name); \
255 geo_base_rep<T,M>::loaded_map().insert (make_pair(name(), base::get_count()));\
257#define _RHEOLEF_set_serial_number(M) \
260geo_basic<T,M>::set_serial_number (size_type i) \
262 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
263 check_macro (ptr != 0, "cannot set_serial_number on geo_domains"); \
264 ptr->set_serial_number(i); \
265 geo_base_rep<T,M>::loaded_map().insert (make_pair(name(), base::get_count()));\
267#define _RHEOLEF_reset_order(M) \
270geo_basic<T,M>::reset_order (size_type order) \
272 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
273 check_macro (ptr != 0, "cannot reset_order on geo_domains"); \
274 ptr->reset_order(order); \
276#define _RHEOLEF_set_nodes(M) \
279geo_basic<T,M>::set_nodes (const disarray<node_type,M>& x) \
281 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
282 check_macro (ptr != 0, "cannot set_nodes on geo_domains"); \
285#define _RHEOLEF_set_coordinate_system(M) \
288geo_basic<T,M>::set_coordinate_system (coordinate_type sys_coord) \
290 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
291 check_macro (ptr != 0, "cannot set_coordinate_system on geo_domains"); \
292 ptr->set_coordinate_system(sys_coord); \
294#define _RHEOLEF_set_dimension(M) \
297geo_basic<T,M>::set_dimension (size_type dim) \
299 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
300 check_macro (ptr != 0, "cannot set_dimension on geo_domains"); \
301 ptr->set_dimension(dim); \
303#define _RHEOLEF_build_by_subdividing(M) \
306geo_basic<T,M>::build_by_subdividing ( \
307 const geo_basic<T,M>& omega, \
310 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
311 check_macro (ptr != 0, "cannot build_by_subdividing on geo_domains"); \
312 ptr->build_by_subdividing (omega, k); \
314#define _RHEOLEF_build_from_data(M) \
317geo_basic<T,M>::build_from_data ( \
318 const geo_header& hdr, \
319 const disarray<node_type, M>& node, \
320 std::array<disarray<geo_element_auto<>,M>, reference_element::max_variant>& tmp_geo_element, \
323 geo_rep<T,M>* ptr = dynamic_cast<geo_rep<T,M>*>(base::pointer()); \
324 check_macro (ptr != 0, "cannot build_from_data on geo_domains"); \
325 ptr->build_from_data (hdr, node, tmp_geo_element, do_upgrade); \
337#ifdef _RHEOLEF_HAVE_MPI
350#undef _RHEOLEF_set_nodes
351#undef _RHEOLEF_reset_order
352#undef _RHEOLEF_set_coordinate_system
353#undef _RHEOLEF_set_dimension
354#undef _RHEOLEF_set_name
355#undef _RHEOLEF_build_from_data
356#undef _RHEOLEF_build_by_subdividing
362#define _RHEOLEF_reset_order(M) \
365geo_rep<T,M>::reset_order (size_type new_order) \
367 if (new_order == base::_piola_basis.degree()) return; \
368 base::_piola_basis.reset_family_index (new_order); \
369 size_type dis_nnod = space_numbering::dis_ndof (base::_piola_basis, base::_gs, base::_gs._map_dimension); \
370 size_type nnod = space_numbering::ndof (base::_piola_basis, base::_gs, base::_gs._map_dimension); \
371 base::_gs.node_ownership = distributor (nnod, base::comm(), nnod); \
372 disarray<point_basic<T>, M> new_node (base::_gs.node_ownership); \
373 for (size_type iv = 0, nv = base::_gs.ownership_by_dimension[0].size(); iv < nv; iv++) { \
374 new_node [iv] = base::_node [iv]; \
376 base::_node = new_node; \
377 build_external_entities (); \
380#ifdef _RHEOLEF_HAVE_MPI
383#undef _RHEOLEF_reset_order
387template <
class T,
class M>
391 if (_name ==
"*nogeo*" || _name ==
"unnamed") {
394#ifdef DEBUG_LOADED_TABLE_TRACE
396 dis_warning_macro (
"geo::~geo(\""<<name()<<
"\"): was not present in the preloaded table");
406template <
class T,
class M>
410 if (_serial_number == 0)
return _name;
412 sprintf (buffer,
"%.3d",
int(_serial_number));
413 return _name +
"-" + buffer;
415template <
class T,
class M>
420 iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
422 if (name == dom.name())
return true;
426template <
class T,
class M>
427const domain_indirect_basic<M>&
431 iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
433 if (name == dom.name())
return dom;
435 error_macro (
"undefined domain \""<<name<<
"\" in mesh \"" << _name <<
"\"");
436 return *(_domains.begin());
438template <
class T,
class M>
443 iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
445 if (dom.name() == curr_dom.name()) {
452 _domains.push_back (dom);
454template <
class T,
class M>
461 sz += _geo_element [variant].size();
465template <
class T,
class M>
472 sz += _geo_element [variant].dis_size();
476template <
class T,
class M>
480 std::vector<size_type>& dis_inod)
const
484template <
class T,
class M>
489 for (
size_type j = 0; j < _dimension; j++) {
490 _xmin[j] = std::numeric_limits<T>::max();
491 _xmax[j] = -std::numeric_limits<T>::max();
493 for (
size_type j = _dimension+1; j < 3; j++) {
494 _xmin [j] = _xmax [j] =
T(0);
496 for (
size_type inod = 0, nnod = _node.size(); inod < nnod; inod++) {
498 for (
size_type j = 0 ; j < _dimension; j++) {
499 _xmin[j] = min(x[j], _xmin[j]);
500 _xmax[j] = max(x[j], _xmax[j]);
504#ifdef _RHEOLEF_HAVE_MPI
506 for (
size_type j = 0 ; j < _dimension; j++) {
507 _xmin[j] = mpi::all_reduce (comm(), _xmin[j], mpi::minimum<T>());
508 _xmax[j] = mpi::all_reduce (comm(), _xmax[j], mpi::maximum<T>());
514 _hmin = std::numeric_limits<T>::max();
516 for (
size_t iedg = 0, nedg = geo_element_ownership(1).size(); iedg < nedg; ++iedg) {
520 T hloc =
dist(x0,x1);
521 _hmin = min(_hmin, hloc);
522 _hmax = max(_hmax, hloc);
524#ifdef _RHEOLEF_HAVE_MPI
526 _hmin = mpi::all_reduce (comm(), _hmin, mpi::minimum<T>());
527 _hmax = mpi::all_reduce (comm(), _hmax, mpi::maximum<T>());
531template <
class T,
class M>
538 last_ige += _geo_element [variant].size();
539 if (ige < last_ige)
return _geo_element [variant] [ige - first_ige];
540 first_ige = last_ige;
542 error_macro (
"geo_element index " << ige <<
" out of range [0:"<< last_ige <<
"[");
543 return _geo_element [0][0];
545template <
class T,
class M>
552 last_ige += _geo_element [variant].size();
553 if (ige < last_ige)
return _geo_element [variant] [ige - first_ige];
554 first_ige = last_ige;
556 error_macro (
"geo_element index " << ige <<
" out of range [0:"<< last_ige <<
"[");
557 return _geo_element [0][0];
563template <
class T,
class M>
568 return _gs.dis_inod2dis_iv (dis_inod);
570template <
class T,
class M>
575 return _gs.dis_iv2dis_inod (dis_iv);
580template <
class T,
class M>
585 if (_gs.ownership_by_dimension[dim].is_owned (dis_ige)) {
586 size_type first_dis_ige = _gs.ownership_by_dimension[dim].first_index();
588 return get_geo_element (dim, ige);
592 size_type dis_igev = _gs.dis_ige2dis_igev_by_dimension (dim, dis_ige, variant);
599template <
class T,
class M>
603 if (omega.have_domain_indirect (
"boundary"))
return;
604 omega.neighbour_guard();
605 size_type map_dim = omega.map_dimension();
606 check_macro (map_dim > 0,
"undefined boundary for 0D geometry");
608 std::vector<size_type> isid_list;
609 for (
size_type isid = 0, nsid = omega.size(sid_dim); isid < nsid; isid++) {
610 const geo_element& S = omega.get_geo_element (sid_dim, isid);
611 if (S.
master(1) != std::numeric_limits<size_type>::max())
continue;
612 isid_list.push_back (isid);
614 communicator comm = omega.sizes().ownership_by_dimension[sid_dim].comm();
615 domain_indirect_basic<M> isid_dom (omega,
"boundary", sid_dim, comm, isid_list);
616 omega.insert_domain_indirect (isid_dom);
618template <
class T,
class M>
622 if (omega.have_domain_indirect (
"internal_sides"))
return;
623 omega.neighbour_guard();
624 size_type map_dim = omega.map_dimension();
625 check_macro (map_dim > 0,
"undefined internal_sides for 0D geometry");
627 std::vector<size_type> isid_list;
628 for (
size_type isid = 0, nsid = omega.size(sid_dim); isid < nsid; isid++) {
629 const geo_element& S = omega.get_geo_element (sid_dim, isid);
630 if (S.
master(1) == std::numeric_limits<size_type>::max())
continue;
631 isid_list.push_back (isid);
633 communicator comm = omega.sizes().ownership_by_dimension[sid_dim].comm();
635 isid_dom.set_broken (
true);
636 omega.insert_domain_indirect (isid_dom);
638template <
class T,
class M>
642 if (omega.have_domain_indirect (
"sides"))
return;
643 omega.neighbour_guard();
644 size_type map_dim = omega.map_dimension();
645 check_macro (map_dim > 0,
"undefined sides for 0D geometry");
647 std::vector<size_type> isid_list;
648 for (
size_type isid = 0, nsid = omega.size(sid_dim); isid < nsid; isid++) {
650 isid_list.push_back (isid);
652 communicator comm = omega.sizes().ownership_by_dimension[sid_dim].comm();
653 domain_indirect_basic<M> isid_dom (omega,
"sides", sid_dim, comm, isid_list);
654 isid_dom.set_broken (
true);
655 omega.insert_domain_indirect (isid_dom);
660#define _RHEOLEF_instanciation(T,M) \
661template class geo_base_rep<T,M>; \
662template class geo_rep<T,M>; \
663template class geo_basic<T,M>; \
664template geo_basic<T,M> geo_load (const std::string& name); \
665template geo_basic<T,M> compact (const geo_basic<T,M>&); \
666template void boundary_guard (const geo_basic<T,M>&); \
667template void internal_sides_guard (const geo_basic<T,M>&); \
668template void sides_guard (const geo_basic<T,M>&);
671#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciation(T, M, A)
field::size_type size_type
see the Float page for the full documentation
see the communicator page for the full documentation
the finite element boundary domain
abstract base interface class
base class for M=sequential or distributed meshes representations
size_type dis_size() const
void dis_inod(const geo_element &K, std::vector< size_type > &dis_inod) const
size_type dis_iv2dis_inod(size_type dis_iv) const
void insert_domain_indirect(const domain_indirect_basic< M > &dom) const
const domain_indirect_basic< M > & get_domain_indirect(size_type i) const
base::size_type size_type
std::unordered_map< std::string, void * > loaded_map_t
const_reference get_geo_element(size_type dim, size_type ige) const
base::const_reference const_reference
const_reference dis_get_geo_element(size_type dim, size_type dis_ige) const
size_type dis_inod2dis_iv(size_type dis_inod) const
static loaded_map_t & loaded_map()
base::reference reference
bool have_domain_indirect(const std::string &name) const
generic mesh with rerefence counting
see the geo_element page for the full documentation
size_type master(bool i) const
base::size_type size_type
base::geo_element_map_type geo_element_map_type
sequential mesh representation
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
#define trace_macro(message)
#define error_macro(message)
#define dis_warning_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_set_name(M)
#define _RHEOLEF_build_by_subdividing(M)
#define _RHEOLEF_set_dimension(M)
#define _RHEOLEF_set_coordinate_system(M)
#define _RHEOLEF_build_from_data(M)
#define _RHEOLEF_reset_order(M)
#define _RHEOLEF_set_serial_number(M)
#define _RHEOLEF_set_nodes(M)
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
This file is part of Rheolef.
string delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
void boundary_guard(const geo_basic< T, M > &omega)
geo_basic< T, M > compact(const geo_basic< T, M > &gamma)
string get_basename(const string &name)
get_basename: see the rheostream page for the full documentation
void sides_guard(const geo_basic< T, M > &omega)
geo_basic< T, M > geo_load(const std::string &filename)
sequential mesh with reference counting
T dist(const point_basic< T > &x, const point_basic< T > &y)
void internal_sides_guard(const geo_basic< T, M > &omega)