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)
43 _u.resize (
_V.iu_ownership(), init_value);
44 _b.resize (
_V.ib_ownership(), init_value);
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);
89 check_macro (
dis_scatch(ids, comp_ios_ownership.
comm(),
"\nfield"),
"read field failed");
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();
132 check_macro (
version <= 3,
"unexpected field version " <<
version);
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();
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;
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());
426 #endif // _RHEOLEF_HAVE_MPI
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) {
481 #ifdef _RHEOLEF_HAVE_MPI
483 mpi::broadcast (mpi::communicator(),
value, ie_proc);
485 #endif // _RHEOLEF_HAVE_MPI
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
533 #endif // _RHEOLEF_HAVE_MPI