rheolef  6.5
field_seq_put_vtk.cc
Go to the documentation of this file.
1 #include "rheolef/field.h"
2 #include "rheolef/field_expr_ops.h"
3 #include "rheolef/field_nonlinear_expr_ops.h"
4 #include "rheolef/interpolate.h"
5 #include "rheolef/piola.h"
6 #include "rheolef/rheostream.h"
7 #include "rheolef/iorheo.h"
8 #include "rheolef/field_evaluate.h"
9 #include "rheolef/space_component.h"
10 namespace rheolef {
11 using namespace std;
12 
13 template <class T> odiststream& geo_put_vtk (odiststream& ods, const geo_basic<T,sequential>& omega,
14  const numbering<T,sequential>& my_numb, const array<point_basic<T>,sequential>& my_node, bool append_data);
15 template <class T> odiststream& geo_put_vtk (odiststream& ods, const geo_basic<T,sequential>& omega,
16  size_t my_order, const array<point_basic<T>,sequential>& my_node, bool append_data);
17 
18 template <class T>
19 static
20 field_basic<T,sequential>
22 {
24  check_macro (Vh.get_geo() == uh.get_geo(), "interpolate: incompatible geos");
26  const space_basic<T,sequential>& Uh = uh.get_space();
27  const geo_basic<T,sequential>& omega = uh.get_geo();
28  std::vector<bool> dof_done (Vh.ndof(), false);
29  std::vector<size_type> u_dis_idof;
30  std::vector<size_type> v_dis_idof;
31  std::vector<point_basic<T> > hat_x;
32  basis_basic<T> fem_basis = Uh.get_numbering().get_basis();
33  basis_basic<T> pointset = Vh.get_numbering().get_basis();
34  basis_on_pointset<T> fem_basis_on_geo_nodes (pointset, fem_basis);
35  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
36  const geo_element& K = omega[ie];
37  Vh.get_numbering().get_basis().hat_node (K.variant(), hat_x);
38  Uh.dis_idof (K, u_dis_idof);
39  Vh.dis_idof (K, v_dis_idof);
40  for (size_type v_loc_idof = 0, v_loc_ndof = v_dis_idof.size(); v_loc_idof < v_loc_ndof; v_loc_idof++) {
41  if (dof_done[v_dis_idof[v_loc_idof]]) continue;
42  dof_done[v_dis_idof[v_loc_idof]] = true;
43  vh.dof (v_dis_idof[v_loc_idof]) = field_evaluate (uh, fem_basis_on_geo_nodes, K, u_dis_idof, v_loc_idof);
44  }
45  }
46  return vh;
47 }
48 template <class T>
49 static
50 field_basic<T,sequential>
52 {
53  if (uh.size() == 0) {
54  return interpolate_scalar (Vh, uh);
55  }
58  for (size_type i_comp = 0, n_comp = uh.size(); i_comp < n_comp; i_comp++) {
59  vh [i_comp] = interpolate_scalar (space_basic<T,sequential>(Vh[i_comp]), field_basic<T,sequential>(uh[i_comp]));
60  }
61  return vh;
62 }
63 template <class T>
64 odiststream&
65 put_vtk_scalar_values (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_header)
66 {
68  ostream& vtk = ods.os();
69  size_type degree = uh.get_space().get_numbering().get_basis().degree();
70  vtk << setprecision(numeric_limits<T>::digits10);
71  if (put_header) {
72  std::string data_type = (degree == 0) ? "CELL_DATA" : "POINT_DATA";
73  vtk << data_type << " " << uh.ndof() << endl;
74  }
75  vtk << "SCALARS " << name << " float" << endl
76  << "LOOKUP_TABLE default" << endl;
77  for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++)
78  vtk << uh.dof(idof) << endl;
79  vtk << endl;
80  return ods;
81 }
82 template <class T>
83 odiststream&
84 put_vtk_vector_values (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_header)
85 {
87  ostream& vtk = ods.os();
88  size_type n_comp = uh.size();
89  std::vector<field_component_const<T,sequential> > uh_comp (n_comp);
90  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
91  uh_comp[i_comp] = uh[i_comp];
92  }
94  field_basic<T,sequential> norm_uh = interpolate (Xh, norm(uh));
95  put_vtk_scalar_values (ods, norm_uh, name+"_norm", put_header);
96  vtk << setprecision(numeric_limits<T>::digits10)
97  << "VECTORS " << name << " float" << endl;
98  std::vector<T> u_dof (n_comp);
99  for (size_type idof = 0, ndof = uh_comp[0].ndof(); idof < ndof; idof++) {
100  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
101  vtk << uh_comp[i_comp].dof (idof);
102  if (i_comp != 2) vtk << " ";
103  }
104  for (size_type i_comp = n_comp; i_comp < 3; i_comp++) {
105  vtk << "0";
106  if (i_comp != 2) vtk << " ";
107  }
108  vtk << endl;
109  }
110  vtk << endl;
111  return ods;
112 }
113 template <class T>
114 odiststream&
115 put_vtk_tensor_values (odiststream& ods, const field_basic<T,sequential>& tau_h, std::string name, bool put_header)
116 {
118  ostream& vtk = ods.os();
119  space_basic<T,sequential> Xh (tau_h.get_geo(), tau_h.get_approx());
120  field_basic<T,sequential> norm_tau_h = interpolate (Xh, norm(tau_h));
121  put_vtk_scalar_values (ods, norm_tau_h, name+"_norm", put_header);
122  vtk << setprecision(numeric_limits<T>::digits10)
123  << "TENSORS " << name << " float" << endl;
124  size_type d = tau_h.get_geo().dimension();
125  switch (d) {
126  case 1: {
127  field_basic<T,sequential> t00 = tau_h(0,0);
128  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
129  vtk << t00.dof(idof) << " 0 0" << endl
130  << "0 0 0" << endl
131  << "0 0 0" << endl;
132  }
133  break;
134  }
135  case 2: {
136  field_basic<T,sequential> t00 = tau_h(0,0);
137  field_basic<T,sequential> t01 = tau_h(0,1);
138  field_basic<T,sequential> t11 = tau_h(1,1);
139  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
140  vtk << t00.dof(idof) << " " << t01.dof(idof) << " 0" << endl
141  << t01.dof(idof) << " " << t11.dof(idof) << " 0" << endl
142  << "0 0 0" << endl;
143  }
144  break;
145  }
146  default: {
147  field_basic<T,sequential> t00 = tau_h(0,0);
148  field_basic<T,sequential> t01 = tau_h(0,1);
149  field_basic<T,sequential> t11 = tau_h(1,1);
150  field_basic<T,sequential> t02 = tau_h(0,2);
151  field_basic<T,sequential> t12 = tau_h(1,2);
152  field_basic<T,sequential> t22 = tau_h(2,2);
153  for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
154  vtk << t00.dof(idof) << " " << t01.dof(idof) << " " << t02.dof(idof) << endl
155  << t01.dof(idof) << " " << t11.dof(idof) << " " << t12.dof(idof) << endl
156  << t02.dof(idof) << " " << t12.dof(idof) << " " << t22.dof(idof) << endl;
157  }
158  }
159  }
160  return ods;
161 }
162 template <class T>
163 odiststream&
164 field_put_vtk (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_geo)
165 {
167  size_type degree = uh.get_space().get_numbering().get_basis().degree();
168  size_type order = uh.get_space().get_geo().order();
170  if (degree == 0) {
171  vh = uh;
172  if (put_geo) {
173  geo_put_vtk (ods, vh.get_geo(), order, vh.get_geo().get_nodes(), false);
174  }
175  } else if (order <= degree && ! uh.get_space().get_numbering().has_compact_support_inside_element()) {
176  vh = uh;
177  if (put_geo) {
178  geo_put_vtk (ods, vh.get_geo(), degree, vh.get_space().get_xdofs(), false);
179  }
180  } else { // order > degree or discontinuous
181  std::string approx = "P" + itos(order);
182  if (uh.get_space().get_numbering().has_compact_support_inside_element()) approx += "d";
183  space_basic<T,sequential> Vh (uh.get_geo(), approx, uh.get_space().valued());
184  vh = interpolate (Vh, uh);
185  if (put_geo) {
186  geo_put_vtk (ods, vh.get_geo(), vh.get_space().get_numbering(), vh.get_space().get_xdofs(), false);
187  }
188  }
189  bool put_header = put_geo;
190  if (name == "") { name = vh.get_space().valued(); }
191  switch (vh.get_space().valued_tag()) {
192  case space_constant::scalar: put_vtk_scalar_values (ods, vh, name, put_header); break;
193  case space_constant::vector: put_vtk_vector_values (ods, vh, name, put_header); break;
194  case space_constant::tensor: put_vtk_tensor_values (ods, vh, name, put_header); break;
195  default: error_macro ("put_vtk: do not known how to print " << vh.valued() << "-valued field");
196  }
197  return ods;
198 }
199 template <class T>
200 odiststream&
202 {
203  return field_put_vtk (ods, uh, "", true);
204 }
205 template odiststream& field_put_vtk<Float> (odiststream&, const field_basic<Float,sequential>&, std::string, bool);
206 template odiststream& field_put_vtk<Float> (odiststream&, const field_basic<Float,sequential>&);
207 
208 }// namespace rheolef
209