rheolef  7.0
field_seq_put_vtk.cc
Go to the documentation of this file.
1 #include "rheolef/field.h"
2 #include "rheolef/piola.h"
3 #include "rheolef/rheostream.h"
4 #include "rheolef/iorheo.h"
5 #include "rheolef/field_evaluate.h"
6 #include "rheolef/space_component.h"
7 #include "rheolef/interpolate.h"
8 
9 namespace rheolef {
10 using namespace std;
11 
12 template <class T> odiststream& geo_put_vtk (odiststream& ods, const geo_basic<T,sequential>& omega,
13  const numbering<T,sequential>& my_numb, const disarray<point_basic<T>,sequential>& my_node, bool append_data);
14 template <class T> odiststream& geo_put_vtk (odiststream& ods, const geo_basic<T,sequential>& omega,
15  size_t my_order, const disarray<point_basic<T>,sequential>& my_node, bool append_data);
16 
17 template <class T>
18 static
19 field_basic<T,sequential>
21 {
23  check_macro (Vh.get_geo() == uh.get_geo(), "interpolate: incompatible geos");
25  const space_basic<T,sequential>& Uh = uh.get_space();
26  const geo_basic<T,sequential>& omega = uh.get_geo();
27  std::vector<bool> dof_done (Vh.ndof(), false);
28  std::vector<size_type> u_dis_idof;
29  std::vector<size_type> v_dis_idof;
30  basis_basic<T> fem_basis = Uh.get_numbering().get_basis();
31  basis_basic<T> pointset = Vh.get_numbering().get_basis();
32  basis_on_pointset<T> fem_basis_on_geo_nodes (pointset, fem_basis);
33  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
34  const geo_element& K = omega[ie];
35  const std::vector<point_basic<T> >& hat_x = Vh.get_numbering().get_basis().hat_node (K.variant());
36  Uh.dis_idof (K, u_dis_idof);
37  Vh.dis_idof (K, v_dis_idof);
38  for (size_type v_loc_idof = 0, v_loc_ndof = v_dis_idof.size(); v_loc_idof < v_loc_ndof; v_loc_idof++) {
39  if (dof_done[v_dis_idof[v_loc_idof]]) continue;
40  dof_done[v_dis_idof[v_loc_idof]] = true;
41  vh.dof (v_dis_idof[v_loc_idof]) = field_evaluate (uh, fem_basis_on_geo_nodes, K, u_dis_idof, v_loc_idof);
42  }
43  }
44  return vh;
45 }
46 template <class T>
47 static
50 {
51  if (uh.size() == 0) {
52  return interpolate_scalar (Vh, uh);
53  }
56  for (size_type i_comp = 0, n_comp = uh.size(); i_comp < n_comp; i_comp++) {
57  vh [i_comp] = interpolate_scalar (space_basic<T,sequential>(Vh[i_comp]), field_basic<T,sequential>(uh[i_comp]));
58  }
59  return vh;
60 }
61 template <class T>
63 put_vtk_scalar_values (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_header)
64 {
66  ostream& vtk = ods.os();
67  size_type degree = uh.get_space().get_numbering().get_basis().degree();
68  vtk << setprecision(numeric_limits<T>::digits10);
69  if (put_header) {
70  std::string data_type = (degree == 0) ? "CELL_DATA" : "POINT_DATA";
71  vtk << data_type << " " << uh.ndof() << endl;
72  }
73  vtk << "SCALARS " << name << " float" << endl
74  << "LOOKUP_TABLE default" << endl;
75  for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++)
76  vtk << uh.dof(idof) << endl;
77  vtk << endl;
78  return ods;
79 }
80 template <class T>
82 put_vtk_vector_values (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_header)
83 {
85  ostream& vtk = ods.os();
86  size_type n_comp = uh.size();
87  std::vector<field_component_const<T,sequential> > uh_comp (n_comp);
88  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
89  uh_comp[i_comp].proxy_assign (uh[i_comp]);
90  }
92  field_basic<T,sequential> norm_uh = interpolate (Xh, norm(uh));
93  put_vtk_scalar_values (ods, norm_uh, name+"_norm", put_header);
94  vtk << setprecision(numeric_limits<T>::digits10)
95  << "VECTORS " << name << " float" << endl;
96  std::vector<T> u_dof (n_comp);
97  for (size_type idof = 0, ndof = uh_comp[0].ndof(); idof < ndof; idof++) {
98  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
99  vtk << uh_comp[i_comp].dof (idof);
100  if (i_comp != 2) vtk << " ";
101  }
102  for (size_type i_comp = n_comp; i_comp < 3; i_comp++) {
103  vtk << "0";
104  if (i_comp != 2) vtk << " ";
105  }
106  vtk << endl;
107  }
108  vtk << endl;
109  return ods;
110 }
111 template <class T>
113 put_vtk_tensor_values (odiststream& ods, const field_basic<T,sequential>& tau_h, std::string name, bool put_header)
114 {
116  ostream& vtk = ods.os();
117  space_basic<T,sequential> Xh (tau_h.get_geo(), tau_h.get_approx());
118  field_basic<T,sequential> norm_tau_h = interpolate (Xh, norm(tau_h));
119  put_vtk_scalar_values (ods, norm_tau_h, name+"_norm", put_header);
120  vtk << setprecision(numeric_limits<T>::digits10)
121  << "TENSORS " << name << " float" << endl;
122  size_type d = tau_h.get_geo().dimension();
123  switch (d) {
124  case 1: {
125  field_basic<T,sequential> t00 = tau_h(0,0);
126  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
127  vtk << t00.dof(idof) << " 0 0" << endl
128  << "0 0 0" << endl
129  << "0 0 0" << endl;
130  }
131  break;
132  }
133  case 2: {
134  field_basic<T,sequential> t00 = tau_h(0,0);
135  field_basic<T,sequential> t01 = tau_h(0,1);
136  field_basic<T,sequential> t11 = tau_h(1,1);
137  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
138  vtk << t00.dof(idof) << " " << t01.dof(idof) << " 0" << endl
139  << t01.dof(idof) << " " << t11.dof(idof) << " 0" << endl
140  << "0 0 0" << endl;
141  }
142  break;
143  }
144  default: {
145  field_basic<T,sequential> t00 = tau_h(0,0);
146  field_basic<T,sequential> t01 = tau_h(0,1);
147  field_basic<T,sequential> t11 = tau_h(1,1);
148  field_basic<T,sequential> t02 = tau_h(0,2);
149  field_basic<T,sequential> t12 = tau_h(1,2);
150  field_basic<T,sequential> t22 = tau_h(2,2);
151  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
152  vtk << t00.dof(idof) << " " << t01.dof(idof) << " " << t02.dof(idof) << endl
153  << t01.dof(idof) << " " << t11.dof(idof) << " " << t12.dof(idof) << endl
154  << t02.dof(idof) << " " << t12.dof(idof) << " " << t22.dof(idof) << endl;
155  }
156  }
157  }
158  return ods;
159 }
160 template <class T>
162 field_put_vtk (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_geo)
163 {
165  size_type degree = uh.get_space().get_numbering().get_basis().degree();
166  size_type order = uh.get_space().get_geo().order();
168  if (degree == 0) {
169  vh = uh;
170  if (put_geo) {
171  geo_put_vtk (ods, vh.get_geo(), order, vh.get_geo().get_nodes(), false);
172  }
173  } else if (order <= degree && ! uh.get_space().get_numbering().has_compact_support_inside_element()) {
174  vh = uh;
175  if (put_geo) {
176  geo_put_vtk (ods, vh.get_geo(), degree, vh.get_space().get_xdofs(), false);
177  }
178  } else { // order > degree or discontinuous
179  std::string approx = "P" + itos(order);
180  if (uh.get_space().get_numbering().has_compact_support_inside_element()) approx += "d";
181  space_basic<T,sequential> Vh (uh.get_geo(), approx, uh.get_space().valued());
182  vh = interpolate (Vh, uh);
183  if (put_geo) {
184  geo_put_vtk (ods, vh.get_geo(), vh.get_space().get_numbering(), vh.get_space().get_xdofs(), false);
185  }
186  }
187  bool put_header = put_geo;
188  if (name == "") { name = vh.get_space().valued(); }
189  switch (vh.get_space().valued_tag()) {
190  case space_constant::scalar: put_vtk_scalar_values (ods, vh, name, put_header); break;
191  case space_constant::vector: put_vtk_vector_values (ods, vh, name, put_header); break;
192  case space_constant::tensor: put_vtk_tensor_values (ods, vh, name, put_header); break;
193  default: error_macro ("put_vtk: do not known how to print " << vh.valued() << "-valued field");
194  }
195  return ods;
196 }
197 template <class T>
200 {
201  return field_put_vtk (ods, uh, "", true);
202 }
203 template odiststream& field_put_vtk<Float> (odiststream&, const field_basic<Float,sequential>&, std::string, bool);
204 template odiststream& field_put_vtk<Float> (odiststream&, const field_basic<Float,sequential>&);
205 
206 }// namespace rheolef
size_type ndof() const
Definition: field.h:318
#define error_macro(message)
Definition: compiler.h:100
void dis_idof(const geo_element &K, std::vector< size_type > &dis_idof) const
Definition: space.h:452
std::size_t size_type
Definition: field.h:234
const basis_basic< T > & get_basis() const
Definition: numbering.h:77
field - piecewise polynomial finite element field
odiststream & geo_put_vtk(odiststream &ods, const geo_basic< T, sequential > &omega, const numbering< T, sequential > &my_numb, const disarray< point_basic< T >, sequential > &my_node, bool append_data)
const numbering< T, sequential > & get_numbering() const
Definition: space.h:382
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
size_type size() const
Definition: field.h:309
std::string get_approx() const
Definition: field.h:286
odiststream & put_vtk_tensor_values(odiststream &ods, const field_basic< T, sequential > &tau_h, std::string name, bool put_header)
STL namespace.
irheostream, orheostream - large data streams
Definition: compiler.h:7
T & dof(size_type idof)
Definition: field.h:553
static field_basic< T, sequential > interpolate(const space_basic< T, sequential > &Vh, const field_basic< T, sequential > &uh)
const std::string & valued() const
Definition: field.h:288
const geo_type & get_geo() const
Definition: field.h:284
odiststream & put_vtk_scalar_values(odiststream &ods, const field_basic< T, sequential > &uh, std::string name, bool put_header)
T norm(const vec< T, M > &x)
Definition: vec.h:341
#define check_macro(ok_condition, message)
Definition: compiler.h:104
size_t d
const geo_basic< T, sequential > & get_geo() const
Definition: space.h:375
void put_header(odiststream &out, const branch_basic< T, sequential > &b)
Definition: branch.cc:30
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
const space_type & get_space() const
Definition: field.h:283
geo_element - element of a mesh
odiststream & put_vtk_vector_values(odiststream &ods, const field_basic< T, sequential > &uh, std::string name, bool put_header)
std::string itos(std::string::size_type i)
static field_basic< T, sequential > interpolate_scalar(const space_basic< T, sequential > &Vh, const field_basic< T, sequential > &uh)
variant_type variant() const
template odiststream & field_put_vtk< Float >(odiststream &, const field_basic< Float, sequential > &, std::string, bool)
std::ostream & os()
Definition: diststream.h:166
size_type size(size_type dim) const
Definition: geo.h:1063
T field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, size_t q)
verbose clean transpose logscale grid shrink ball stereo iso volume velocity elevation bezieradapt label color color format format format format format format format format format format format format format format format format vtk