21#include "rheolef/field_wdof_sliced.h"
22#include "rheolef/field.h"
23#include "rheolef/rheostream.h"
24#include "rheolef/iorheo.h"
25#include "rheolef/piola_util.h"
26#include "rheolef/field_expr.h"
34template <
class T,
class M>
41 _dis_dof_indexes_requires_update(true),
42 _dis_dof_assembly_requires_update(true)
47template <
class T,
class M>
55 _u.resize (_V.iu_ownership(), init_value);
56 _b.resize (_V.ib_ownership(), init_value);
57 dis_dof_indexes_requires_update();
58 dis_dof_assembly_requires_update();
63template <
class T,
class M>
72 bool read_header =
false)
82 const distributor& comp_ios_ownership = terminal_constit.ios_ownership();
87 distributor comp_ios_ownership (comp_dis_ndof, ids.
comm(), comp_ios_ndof);
92 std::string geo_name, approx;
93 ids >> version >> dis_size1 >> geo_name >> approx;
96 vec<T,M> u_comp_io (comp_ios_ownership, std::numeric_limits<T>::max());
97 u_comp_io.get_values (ids);
99 for (
size_type comp_ios_idof = 0; comp_ios_idof < comp_ios_ndof; comp_ios_idof++) {
100 size_type ios_idof = comp_start_ios_idof + comp_ios_idof;
101 assert_macro (ios_idof < ios_ndof,
"ios_idof="<<ios_idof<<
" out of range [0:"<<ios_ndof<<
"[");
102 u_io [ios_idof] = u_comp_io [comp_ios_idof];
104 comp_start_dis_idof += comp_dis_ndof;
105 comp_start_ios_idof += comp_ios_ndof;
109 typedef typename space_constitution<T,M>::hierarchy_type hier_t;
111 for (
typename hier_t::const_iterator iter = hier_constit.begin(), last = hier_constit.end(); iter != last; ++iter) {
112 const space_constitution<T,M>& curr_constit = *iter;
113 get_field_recursive (ids, u_io, curr_constit, comp_start_dis_idof, comp_start_ios_idof,
true);
117template<
class T,
class M>
120template <
class T,
class M>
125 dis_dof_indexes_requires_update();
133 check_macro (version <= 3,
"unexpected field version " << version);
135 std::string constit_name_input;
139 std::string geo_name, approx;
144 if (_V.get_geo().name() == geo_name) {
145 omega = _V.get_geo();
155 bool have_constit =
false;
158 check_macro (label ==
"header",
"field file format version "<< version <<
": \"header\" keyword not found");
161 if (label ==
"end") {
163 }
else if (label ==
"size") {
166 }
else if (label ==
"constitution") {
167 ids >> constit_name_input;
168 std::istringstream istrstr (constit_name_input);
171 idiststrstr >> constit;
177 error_macro (
"unexpected field header member: \""<<label<<
"\"");
181 check_macro (label ==
"header",
"field file format version "<< version <<
": \"end header\" keyword not found");
182 check_macro (have_constit,
"field file format version "<< version <<
": \"constitution\" keyword not found");
185 if (!(_V.get_constitution() == constit)) {
189 size_type dis_ndof = _V.ownership().dis_size();
192 vec<T,M> u_io (_V.ios_ownership(), std::numeric_limits<T>::max());
193 vec<T,M> u_dof (_V.ownership(), std::numeric_limits<T>::max());
196 u_io.get_values (ids);
197 for (
size_type ios_idof = 0, ios_ndof = _V.ios_ownership().size(); ios_idof < ios_ndof; ios_idof++) {
198 const T& value = u_io [ios_idof];
199 size_type dis_idof = _V.ios_idof2dis_idof (ios_idof);
200 u_dof.dis_entry (dis_idof) = value;
203 u_dof.dis_entry_assembly();
205 bool need_old2new_convert
209 if (!need_old2new_convert) {
210 for (
size_type idof = 0, ndof = _V.ownership().size(); idof < ndof; idof++) {
211 dof(idof) = u_dof [idof];
216 std::string valued = constit.
valued();
217 std::string approx = constit[0].
get_basis().name();
218 std::string geo_name = constit.
get_geo().name();
220 std::string new_constit_name_input = valued +
"(" + approx +
"){" + geo_name +
"}";
221 std::istringstream new_istrstr (new_constit_name_input);
224 new_idiststrstr >> new_constit;
229 for (
size_type i_comp = 0; i_comp < n_comp; i_comp++) {
230 for (
size_type comp_idof = 0; comp_idof < comp_ndof; comp_idof++) {
231 size_type new_idof = comp_idof*n_comp + i_comp;
232 size_type old_idof = i_comp*comp_ndof + comp_idof;
233 dof(new_idof) = u_dof [old_idof];
241template <
class T,
class M>
250 bool write_header =
false)
263 distributor comp_ios_ownership (comp_dis_ndof, comm, (my_proc == io_proc ? comp_dis_ndof : 0));
264 vec<T,M> comp_u_io (comp_ios_ownership, std::numeric_limits<T>::max());
265 for (
size_type comp_idof = 0; comp_idof < comp_ndof; comp_idof++) {
266 size_type idof = comp_start_idof + comp_idof;
267 T value = uh.
dof (idof);
269 assert_macro (ios_dis_idof >= comp_start_dis_idof,
"invalid comp ios index");
270 size_type comp_ios_dis_idof = ios_dis_idof - comp_start_dis_idof;
271 comp_u_io.dis_entry (comp_ios_dis_idof) = value;
273 comp_u_io.dis_entry_assembly();
276 ods << setprecision(std::numeric_limits<Float>::digits10);
278 ods <<
"field" << endl
279 <<
"1 " << comp_dis_ndof << endl
280 << terminal_constit.
get_geo().name() << endl
281 << terminal_constit.
get_basis().name() << endl
285 << setprecision(old_prec);
286 comp_start_idof += comp_ndof;
287 comp_start_dis_idof += comp_dis_ndof;
288 if (comp_start_dis_idof != uh.
dis_ndof()) {
296 for (
typename hier_t::const_iterator iter = hier_constit.begin(), last = hier_constit.end(); iter != last; ++iter) {
298 put_field_recursive (ods, uh, curr_constit, comp_start_idof, comp_start_dis_idof,
false);
301template <
class T,
class M>
306 bool need_header = (get_space().get_constitution().is_hierarchical());
309 ods <<
"field" << endl
312 <<
" constitution " << get_space().get_constitution().name() << endl
313 <<
" size " << get_space().get_constitution().dis_ndof() << endl
314 <<
"end header" << endl
319 ods <<
"field" << endl
320 <<
"1 " << dis_ndof() << endl
322 << terminal_constit.
get_basis().name() << endl
327 put_field_recursive (ods, *
this, get_space().get_constitution(), comp_start_idof, comp_start_dis_idof);
343template <
class T,
class M>
351struct field_put<
T,sequential> {
352 odiststream& operator() (odiststream& ods,
const field_basic<T,sequential>& uh)
const
361 return uh.put_field (ods);
365template <
class T,
class M>
369 field_put<T,M> put_fct;
370 return put_fct (ods, *
this);
375template <
class T,
class M>
380 if (_dis_dof_indexes_requires_update || _dis_dof_assembly_requires_update) {
382 check_macro (nproc == 1,
"field::dis_dof_update() need to be called before field::dis_dof(dis_idof)");
385 if (ownership().is_owned (dis_idof)) {
386 size_type first_idof = ownership().first_index ();
391 space_pair_type blk_dis_iub = _V.data()._idof2blk_dis_iub.dis_at (dis_idof);
393 const T& value = (! blk_dis_iub.
is_blocked()) ? _u.dis_at (dis_iub) : _b.dis_at (dis_iub);
396template <
class T,
class M>
400 dis_dof_assembly_requires_update();
401 const space_pair_type& blk_dis_iub = _V.data()._idof2blk_dis_iub.dis_at (dis_idof);
404 return _u.dis_entry (dis_iub);
406 return _b.dis_entry (dis_iub);
412template <
class T,
class M>
418 std::vector<size_type> dis_idof1 (loc_ndof);
419 _V.dis_idof (K, dis_idof1);
421 std::vector<T> dof (loc_ndof);
422 for (
size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
423 dof [loc_idof] = dis_dof (dis_idof1 [loc_idof]);
428 Eigen::Matrix<T,Eigen::Dynamic,1> b_value (loc_ndof);
429 b.evaluate (K, hat_x, b_value);
432 for (
size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
433 value += dof [loc_idof] * b_value[loc_idof];
437template <
class T,
class M>
444 check_macro (dis_ie != std::numeric_limits<size_type>::max(),
"x="<<x<<
" is outside the domain");
446 T value = std::numeric_limits<T>::max();
447 size_type map_dim = omega.map_dimension();
448 const distributor& ownership = omega.geo_element_ownership(map_dim);
451 std::vector<size_type> dis_inod;
452 if (my_proc == ie_proc) {
456 omega.dis_inod (K, dis_inod);
458 value =
evaluate (K, hat_x, i_comp);
460#ifdef _RHEOLEF_HAVE_MPI
462 mpi::broadcast (mpi::communicator(), value, ie_proc);
467template <
class T,
class M>
478template <
class T,
class M>
486template <
class T,
class M>
487details::field_wdof_sliced<field_basic<T,M>>
492 return details::field_wdof_sliced<field_basic<T,M>> (*
this, ij_comp);
498#define _RHEOLEF_instanciation_base(T,M) \
499template class field_basic<T,M>; \
500template odiststream& operator<< (odiststream&, const field_basic<T,M>&);
503#define _RHEOLEF_instanciation(T,M) \
504_RHEOLEF_instanciation_base(T,M) \
505_RHEOLEF_instanciation_base(std::complex<T>,M)
507#define _RHEOLEF_instanciation(T,M) \
508_RHEOLEF_instanciation_base(T,M)
512#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
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 T & dis_dof(size_type dis_idof) const
size_type dis_ndof() const
idiststream & get(idiststream &ips)
point_basic< T > dis_vector_evaluate(const point_basic< T > &x) const
odiststream & put(odiststream &ops) const
dis_reference dis_dof_entry(size_type dis_idof)
odiststream & put_field(odiststream &ops) const
typename vec< scalar_type, memory_type >::dis_reference dis_reference
T operator()(const point_basic< T > &x) const
const space_type & get_space() const
T evaluate(const geo_element &K, const point_basic< T > &hat_xq, size_type i_comp=0) const
void resize(const space_type &V, const T &init_value=std::numeric_limits< T >::max())
T dis_evaluate(const point_basic< T > &x, size_type i_comp=0) const
see the geo_element page for the full documentation
variant_type variant() const
idiststream: see the diststream page for the full documentation
static size_type io_proc()
This routine returns the rank of a process that can perform i/o.
const communicator & comm() const
std::bitset< last > flag_type
static flag_type format_field
odiststream: see the diststream page for the full documentation
static size_type io_proc()
const basis_basic< T > & get_basis() const
const geo_basic< T, M > & get_geo() const
const basis_basic< T > & get_basis() const
size_type dis_ndof() const
const valued_type & valued_tag() const
rep::hierarchy_type hierarchy_type
const space_constitution_terminal< T, M > & get_terminal() const
const hierarchy_type & get_hierarchy() const
size_type ios_ndof() const
communicator comm() const
bool is_hierarchical() const
const std::string & valued() const
const geo_basic< T, M > & get_geo() const
see the vec page for the full documentation
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)")
size_type tensor_index(valued_type valued_tag, coordinate_type sys_coord, size_type i, size_type j)
This file is part of Rheolef.
void space_constitution_old_get(idiststream &ids, space_constitution< T, M > &constit)
odiststream & visu_vtk_paraview(odiststream &, const field_basic< T, sequential > &)
odiststream & field_put_gmsh_pos(odiststream &, const field_basic< T, sequential > &)
odiststream & visu_gmsh(odiststream &, const field_basic< T, sequential > &)
bool dis_scatch(idiststream &ips, const communicator &comm, std::string ch)
distributed version of scatch(istream&,string)
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)
odiststream & field_put_bamg_bb(odiststream &, const field_basic< T, sequential > &)
odiststream & visu_gnuplot(odiststream &, const field_basic< T, sequential > &)
odiststream & field_put_gmsh(odiststream &, const field_basic< T, sequential > &, std::string)
space_constant::valued_type valued_tag() const
size_type dis_iub() const