1 #include "rheolef/field.h"
22 #include "rheolef/rheostream.h"
23 #include "rheolef/iorheo.h"
24 #include "rheolef/piola_util.h"
25 #include "rheolef/field_expr.h"
33 template <
class T,
class M>
40 _dis_dof_indexes_requires_update(true),
41 _dis_dof_assembly_requires_update(true)
46 template <
class T,
class M>
54 _u.resize (_V.iu_ownership(), init_value);
55 _b.resize (_V.ib_ownership(), init_value);
56 dis_dof_indexes_requires_update();
57 dis_dof_assembly_requires_update();
62 template <
class T,
class M>
71 bool read_header =
false)
81 const distributor& comp_ios_ownership = scalar_constit.ios_ownership();
86 distributor comp_ios_ownership (comp_dis_ndof, ids.comm(), comp_ios_ndof);
91 std::string geo_name, approx;
92 ids >>
version >> dis_size1 >> geo_name >> approx;
95 vec<T,M> u_comp_io (comp_ios_ownership, std::numeric_limits<T>::max());
96 u_comp_io.get_values (ids);
98 for (
size_type comp_ios_idof = 0; comp_ios_idof < comp_ios_ndof; comp_ios_idof++) {
99 size_type ios_idof = comp_start_ios_idof + comp_ios_idof;
100 assert_macro (ios_idof < ios_ndof,
"ios_idof="<<ios_idof<<
" out of range [0:"<<ios_ndof<<
"[");
101 u_io [ios_idof] = u_comp_io [comp_ios_idof];
103 comp_start_dis_idof += comp_dis_ndof;
104 comp_start_ios_idof += comp_ios_ndof;
110 for (
typename hier_t::const_iterator iter = hier_constit.begin(), last = hier_constit.end(); iter != last; ++iter) {
111 const space_constitution<T,M>& curr_constit = *iter;
112 get_field_recursive (ids, u_io, curr_constit, comp_start_dis_idof, comp_start_ios_idof,
true);
116 template<
class T,
class M>
119 template <
class T,
class M>
124 dis_dof_indexes_requires_update();
125 communicator comm = ids.comm();
134 std::string constit_name_input;
138 std::string geo_name, approx;
143 if (_V.get_geo().name() == geo_name) {
144 omega = _V.get_geo();
154 bool have_constit =
false;
157 check_macro (label ==
"header",
"field file format version "<<
version <<
": \"header\" keyword not found");
160 if (label ==
"end") {
162 }
else if (label ==
"size") {
165 }
else if (label ==
"constitution") {
166 ids >> constit_name_input;
167 std::istringstream istrstr (constit_name_input);
168 idiststream idiststrstr (istrstr);
170 idiststrstr >> constit;
176 error_macro (
"unexpected field header member: \""<<label<<
"\"");
180 check_macro (label ==
"header",
"field file format version "<<
version <<
": \"end header\" keyword not found");
181 check_macro (have_constit,
"field file format version "<<
version <<
": \"constitution\" keyword not found");
184 if (!(_V.get_constitution() == constit)) {
191 vec<T,M> u_io (_V.ios_ownership(), std::numeric_limits<T>::max());
192 vec<T,M> u_dof (_V.ownership(), std::numeric_limits<T>::max());
195 u_io.get_values (ids);
196 for (
size_type ios_idof = 0, ios_ndof = _V.ios_ownership().size(); ios_idof < ios_ndof; ios_idof++) {
197 const T& value = u_io [ios_idof];
202 u_dof.dis_entry_assembly();
204 bool need_old2new_convert
208 if (!need_old2new_convert) {
210 dof(idof) = u_dof [idof];
215 std::string valued = constit.
valued();
216 std::string approx = constit[0].
get_basis().name();
217 std::string geo_name = constit.
get_geo().name();
219 std::string new_constit_name_input = valued +
"(" + approx +
"){" + geo_name +
"}";
220 std::istringstream new_istrstr (new_constit_name_input);
221 idiststream new_idiststrstr (new_istrstr);
223 new_idiststrstr >> new_constit;
228 for (
size_type i_comp = 0; i_comp < n_comp; i_comp++) {
229 for (
size_type comp_idof = 0; comp_idof < comp_ndof; comp_idof++) {
230 size_type new_idof = comp_idof*n_comp + i_comp;
231 size_type old_idof = i_comp*comp_ndof + comp_idof;
232 dof(new_idof) = u_dof [old_idof];
240 template <
class T,
class M>
243 put_field_recursive (
249 bool write_header =
false)
258 communicator comm = constit.
comm();
260 size_type io_proc = odiststream::io_proc();
262 distributor comp_ios_ownership (comp_dis_ndof, comm, (my_proc == io_proc ? comp_dis_ndof : 0));
263 vec<T,M> comp_u_io (comp_ios_ownership, std::numeric_limits<T>::max());
264 for (
size_type comp_idof = 0; comp_idof < comp_ndof; comp_idof++) {
265 size_type idof = comp_start_idof + comp_idof;
266 T value = uh.
dof (idof);
268 assert_macro (ios_dis_idof >= comp_start_dis_idof,
"invalid comp ios index");
269 size_type comp_ios_dis_idof = ios_dis_idof - comp_start_dis_idof;
270 comp_u_io.dis_entry (comp_ios_dis_idof) = value;
272 comp_u_io.dis_entry_assembly();
275 ods << setprecision(std::numeric_limits<Float>::digits10);
277 ods <<
"field" << endl
278 <<
"1 " << comp_dis_ndof << endl
279 << scalar_constit.
get_geo().name() << endl
280 << scalar_constit.
get_basis().name() << endl
284 << setprecision(old_prec);
285 comp_start_idof += comp_ndof;
286 comp_start_dis_idof += comp_dis_ndof;
287 if (comp_start_dis_idof != uh.
dis_ndof()) {
295 for (
typename hier_t::const_iterator iter = hier_constit.begin(), last = hier_constit.end(); iter != last; ++iter) {
296 const space_constitution<T,M>& curr_constit = *iter;
297 put_field_recursive (ods, uh, curr_constit, comp_start_idof, comp_start_dis_idof,
false);
300 template <
class T,
class M>
305 bool need_header = (get_space().get_constitution().is_hierarchical());
308 ods <<
"field" << endl
311 <<
" constitution " << get_space().get_constitution().name() << endl
312 <<
" size " << get_space().get_constitution().dis_ndof() << endl
313 <<
"end header" << endl
318 ods <<
"field" << endl
320 << scalar_constit.
get_geo().name() << endl
321 << scalar_constit.
get_basis().name() << endl
326 put_field_recursive (ods, *
this, get_space().get_constitution(), comp_start_idof, comp_start_dis_idof);
342 template <
class T,
class M>
353 iorheo::flag_type format = iorheo::flags(ods.
os()) & iorheo::format_field;
364 template <
class T,
class M>
368 field_put<T,M> put_fct;
369 return put_fct (ods, *
this);
374 template <
class T,
class M>
379 if (_dis_dof_indexes_requires_update || _dis_dof_assembly_requires_update) {
381 check_macro (nproc == 1,
"field::dis_dof_update() need to be called before field::dis_dof(dis_idof)");
384 if (ownership().is_owned (
dis_idof)) {
385 size_type first_idof = ownership().first_index ();
392 const T& value = (! blk_dis_iub.
is_blocked()) ? _u.dis_at (dis_iub) : _b.dis_at (dis_iub);
395 template <
class T,
class M>
399 dis_dof_assembly_requires_update();
403 return _u.dis_entry (dis_iub);
405 return _b.dis_entry (dis_iub);
408 template <
class T,
class M>
412 #ifdef _RHEOLEF_HAVE_MPI
413 std::size_t nproc = ownership().
comm().size();
415 std::size_t do_assembly = mpi::all_reduce (ownership().comm(),
size_t(_dis_dof_assembly_requires_update), std::plus<std::size_t>());
416 std::size_t do_indexes = mpi::all_reduce (ownership().comm(),
size_t(_dis_dof_indexes_requires_update), std::plus<std::size_t>());
418 _u.dis_entry_assembly();
419 _b.dis_entry_assembly();
422 _u.set_dis_indexes (_V.ext_iu_set());
423 _b.set_dis_indexes (_V.ext_ib_set());
427 _dis_dof_indexes_requires_update =
false;
428 _dis_dof_assembly_requires_update =
false;
433 template <
class T,
class M>
439 std::vector<size_type> dis_idof1 (loc_ndof);
440 _V.dis_idof (K, dis_idof1);
442 std::vector<T> dof (loc_ndof);
443 for (
size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
444 dof [loc_idof] = dis_dof (dis_idof1 [loc_idof]);
449 Eigen::Matrix<T,Eigen::Dynamic,1> b_value (loc_ndof);
450 b.evaluate (K, hat_x, b_value);
453 for (
size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
454 value += dof [loc_idof] * b_value[loc_idof];
458 template <
class T,
class M>
465 check_macro (dis_ie != std::numeric_limits<size_type>::max(),
"x="<<x<<
" is outside the domain");
467 T value = std::numeric_limits<T>::max();
473 if (my_proc == ie_proc) {
479 value =
evaluate (K, hat_x, i_comp);
481 #ifdef _RHEOLEF_HAVE_MPI
483 mpi::broadcast (mpi::communicator(), value, ie_proc);
488 template <
class T,
class M>
498 template <
class T,
class M>
506 template <
class T,
class M>
517 #define _RHEOLEF_instanciation_base(T,M) \
518 template class field_basic<T,M>; \
519 template odiststream& operator<< (odiststream&, const field_basic<T,M>&);
522 #define _RHEOLEF_instanciation(T,M) \
523 _RHEOLEF_instanciation_base(T,M) \
524 _RHEOLEF_instanciation_base(std::complex<T>,M)
526 #define _RHEOLEF_instanciation(T,M) \
527 _RHEOLEF_instanciation_base(T,M)
531 #ifdef _RHEOLEF_HAVE_MPI
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
field::size_type size_type
see the Float page for the full documentation
see the distributor page for the full documentation
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type size(size_type iproc) const
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
std::allocator< int >::size_type size_type
const communicator_type & comm() const
const space_type & get_space() const
size_type dis_ndof() const
odiststream & put_field(odiststream &ops) const
const communicator & comm() const
void resize(const space_type &V, const T &init_value=std::numeric_limits< T >::max())
generic mesh with rerefence counting
see the geo_element page for the full documentation
variant_type variant() const
odiststream: see the diststream page for the full documentation
const space_scalar_constitution< T, M > & get_scalar() const
size_type dis_ndof() const
const valued_type & valued_tag() const
const geo_basic< T, M > & get_geo() const
rep::hierarchy_type hierarchy_type
const std::string & valued() const
size_type ios_ndof() const
communicator comm() const
bool is_hierarchical() const
const hierarchy_type & get_hierarchy() const
const basis_basic< T > & get_basis() const
const geo_basic< T, M > & get_geo() const
const basis_basic< T > & get_basis() const
void resize(const distributor &ownership, const T &init_val=std::numeric_limits< T >::max())
#define assert_macro(ok_condition, message)
#define error_macro(message)
#define fatal_macro(message)
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format bamg
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format format format format paraview
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format format gnuplot
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format gmsh_pos
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format gmsh
#define _RHEOLEF_instanciation(T, M)
size_type tensor_index(valued_type valued_tag, coordinate_type sys_coord, size_type i, size_type j)
valued_type valued_tag(const std::string &name)
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)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
size_type dis_ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
odiststream & field_put_gmsh_pos(odiststream &, const field_basic< T, sequential > &)
odiststream & field_put_gmsh(odiststream &, const field_basic< T, sequential > &, std::string)
void space_constitution_old_get(idiststream &ids, space_constitution< T, M > &constit)
odiststream & field_put_bamg_bb(odiststream &, const field_basic< T, sequential > &)
void evaluate(const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, 1 > &value) const
point_basic< T > inverse_piola_transformation(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod, const point_basic< T > &x)
bool dis_scatch(idiststream &ips, const communicator &comm, std::string ch)
distributed version of scatch(istream&,string)
odiststream & visu_gnuplot(odiststream &, const field_basic< T, sequential > &)
odiststream & visu_vtk_paraview(odiststream &, const field_basic< T, sequential > &)
odiststream & visu_gmsh(odiststream &, const field_basic< T, sequential > &)
size_type dis_iub() const