Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_seq_put_vtk.cc
Go to the documentation of this file.
1
21//
22// vtk output
23//
24// author: Pierre.Saramito@imag.fr
25//
26// date: 12 may 1997 update: 23 oct 2011
27//
28#include "rheolef/field.h"
29#include "rheolef/piola_util.h"
30#include "rheolef/rheostream.h"
31#include "rheolef/iorheo.h"
32#include "rheolef/field_evaluate.h"
33#include "rheolef/space_component.h"
34#include "rheolef/diststream.h"
35#include "rheolef/interpolate.h"
36
37#include "geo_seq_put_vtk.h"
38
39namespace rheolef {
40using namespace std;
41
42// ----------------------------------------------------------------------------
43// low order elements : {P0,P1,P1d}-iso-P1
44// when approx = P0 or P1, P1d, paraview is more efficient with it
45// than with the general nonlinear Lagrange elements
46// ----------------------------------------------------------------------------
47template <class T>
48odiststream&
50{
52 ostream& vtk = ods.os();
53 size_type degree = uh.get_space().get_basis().degree();
54 vtk << setprecision(numeric_limits<T>::digits10);
55 if (put_header) {
56 std::string data_type = (degree == 0) ? "CELL_DATA" : "POINT_DATA";
57 vtk << data_type << " " << uh.ndof() << endl;
58 }
59 vtk << "SCALARS " << name << " float" << endl
60 << "LOOKUP_TABLE default" << endl;
61 for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
62 vtk << uh.dof(idof) << endl;
63 }
64 vtk << endl;
65 return ods;
66}
67template <class T>
68odiststream&
70{
72 ostream& vtk = ods.os();
73 // 1. output norm(uh)
74 const basis_basic<T>& b_comp = uh.get_space().get_constitution().get_basis()[0];
75 space_basic<T,sequential> Xh (uh.get_geo(), b_comp.name());
76 field_basic<T,sequential> norm_uh = interpolate (Xh, norm(uh));
77 put_vtk_scalar_values (ods, norm_uh, name+"_norm", put_header);
78 // 1. output uh components
79 vtk << setprecision(numeric_limits<T>::digits10)
80 << "VECTORS " << name << " float" << endl;
81 size_type n_comp = uh.get_space().get_basis().size();
82 for (size_type i_comp_dof = 0, n_comp_dof = uh.ndof()/n_comp; i_comp_dof < n_comp_dof; i_comp_dof++) {
83 for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
84 size_type idof = i_comp_dof*n_comp + i_comp;
85 vtk << uh.dof (idof);
86 if (i_comp != 2) vtk << " ";
87 }
88 for (size_type i_comp = n_comp; i_comp < 3; i_comp++) {
89 vtk << "0";
90 if (i_comp != 2) vtk << " ";
91 }
92 vtk << endl;
93 }
94 vtk << endl;
95 return ods;
96}
97template <class T>
98odiststream&
99put_vtk_tensor_values (odiststream& ods, const field_basic<T,sequential>& tau_h, std::string name, bool put_header)
100{
102 ostream& vtk = ods.os();
103 const basis_basic<T>& b_comp = tau_h.get_space().get_constitution().get_basis()[0];
104 space_basic<T,sequential> Xh (tau_h.get_geo(), b_comp.name());
105 field_basic<T,sequential> norm_tau_h = interpolate (Xh, norm(tau_h));
106 put_vtk_scalar_values (ods, norm_tau_h, name+"_norm", put_header);
107 vtk << setprecision(numeric_limits<T>::digits10)
108 << "TENSORS " << name << " float" << endl;
109 size_type d = tau_h.get_geo().dimension();
110 switch (d) {
111 case 1: {
112 field_basic<T,sequential> t00 = tau_h(0,0);
113 for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
114 vtk << t00.dof(idof) << " 0 0" << endl
115 << "0 0 0" << endl
116 << "0 0 0" << endl;
117 }
118 break;
119 }
120 case 2: {
121 space_constant::coordinate_type sys_coord = tau_h.get_geo().coordinate_system();
122 size_type i_comp00 = space_constant::tensor_index (tau_h.valued_tag(), sys_coord, 0, 0);
123 size_type i_comp01 = space_constant::tensor_index (tau_h.valued_tag(), sys_coord, 0, 1);
124 size_type i_comp11 = space_constant::tensor_index (tau_h.valued_tag(), sys_coord, 1, 1);
125 size_type n_comp = tau_h.get_space().get_basis().size();
126 for (size_type i_comp_dof = 0, n_comp_dof = tau_h.ndof()/n_comp; i_comp_dof < n_comp_dof; i_comp_dof++) {
127 size_type idof00 = i_comp_dof*n_comp + i_comp00;
128 size_type idof01 = i_comp_dof*n_comp + i_comp01;
129 size_type idof11 = i_comp_dof*n_comp + i_comp11;
130 vtk << tau_h.dof(idof00) << " " << tau_h.dof(idof01) << " 0" << endl
131 << tau_h.dof(idof01) << " " << tau_h.dof(idof11) << " 0" << endl
132 << "0 0 0" << endl;
133 }
134 break;
135 }
136 default: {
137 field_basic<T,sequential> t00 = tau_h(0,0);
138 field_basic<T,sequential> t01 = tau_h(0,1);
139 field_basic<T,sequential> t11 = tau_h(1,1);
140 field_basic<T,sequential> t02 = tau_h(0,2);
141 field_basic<T,sequential> t12 = tau_h(1,2);
142 field_basic<T,sequential> t22 = tau_h(2,2);
143 for (size_type idof = 0, ndof = t00.ndof(); idof < ndof; idof++) {
144 vtk << t00.dof(idof) << " " << t01.dof(idof) << " " << t02.dof(idof) << endl
145 << t01.dof(idof) << " " << t11.dof(idof) << " " << t12.dof(idof) << endl
146 << t02.dof(idof) << " " << t12.dof(idof) << " " << t22.dof(idof) << endl;
147 }
148 }
149 }
150 return ods;
151}
152template <class T>
153odiststream&
154field_put_vtk (odiststream& ods, const field_basic<T,sequential>& uh, std::string name, bool put_geo)
155{
157 // 1) output geo
158 basis_basic<T> my_numb = uh.get_space().get_basis();
159 basis_basic<T> my_piola = uh.get_geo().get_piola_basis();
160 size_type degree = my_numb.degree();
161 size_type order = uh.get_space().get_geo().order();
162 size_type subdivide = iorheo::getsubdivide (dout.os());
164 if (subdivide == 0) {
165 if (degree == 0) {
166 // P0
167 vh = uh;
168 if (put_geo) {
169 geo_put_vtk (ods, vh.get_geo(), my_piola, vh.get_geo().get_nodes(), false, vh.get_geo().map_dimension());
170 }
171 } else if (order <= degree && uh.get_space().get_basis().family_name() == "P") {
172 // Pk[d] and order <= degree
173 vh = uh;
174 if (put_geo) {
175 geo_put_vtk (ods, vh.get_geo(), my_numb, vh.get_space().get_xdofs(), false, vh.get_geo().map_dimension());
176 }
177 } else {
178 // not Pk[d] or order > degree
179 // => re-interpolate uh on the finer mesh lattice (eg P1d field on a P2 mesh)
180 size_type k = std::max (order, degree);
181 std::string approx = "P" + std::to_string(k);
182 if (uh.get_space().get_basis().have_compact_support_inside_element()) approx += "d";
183 space_basic<T,sequential> Vh (uh.get_geo(), approx, uh.get_space().valued());
184 warning_macro ("reinterpolate \""<< approx << "\" since \"" << uh.get_approx() << "\" is not Pk[d] or mesh order > degree");
185 vh = interpolate (Vh, uh);
186 if (put_geo) {
187 geo_put_vtk (ods, vh.get_geo(), vh.get_space().get_basis(), vh.get_space().get_xdofs(), false, vh.get_geo().map_dimension());
188 }
189 }
190 } else { // subdivide > 0
192 omega_s.build_by_subdividing (uh.get_geo(), subdivide);
193 space_basic<T,sequential> Vh_s (omega_s, uh.get_space().get_approx());
194 vh = interpolate (Vh_s, uh);
195 // TODO: interpolate: not efficient (octree...)
196 // + when discontinuous, force continuous junctions
197 // => see field_seq_visu_gnuplot.cc for direct output
198 if (put_geo) {
199 geo_put_vtk (ods, vh.get_geo(), vh.get_space().get_basis(), vh.get_space().get_xdofs(), false, vh.get_geo().map_dimension());
200 }
201 }
202 // 2) output values
203 bool put_header = put_geo;
204 if (name == "") { name = vh.get_space().valued(); }
205 switch (vh.get_space().valued_tag()) {
206 case space_constant::scalar: put_vtk_scalar_values (ods, vh, name, put_header); break;
207 case space_constant::vector: put_vtk_vector_values (ods, vh, name, put_header); break;
208 case space_constant::tensor: put_vtk_tensor_values (ods, vh, name, put_header); break;
209 default: error_macro ("put_vtk: do not known how to print " << vh.valued() << "-valued field");
210 }
211 return ods;
212}
213// ----------------------------------------------------------------------------
214// main call
215// ----------------------------------------------------------------------------
216template <class T>
217odiststream&
219{
220 return field_put_vtk (ods, uh, "", true);
221}
222// ----------------------------------------------------------------------------
223// instanciation in library
224// ----------------------------------------------------------------------------
227
228}// namespace rheolef
field::size_type size_type
Definition branch.cc:430
std::string name() const
Definition basis.h:721
size_type degree() const
Definition basis.h:728
valued_type valued_tag() const
Definition field.h:273
std::string get_approx() const
Definition field.h:272
T & dof(size_type idof)
Definition field.h:738
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
size_type ndof() const
Definition field.h:298
const std::string & valued() const
Definition field.h:274
std::size_t size_type
Definition field.h:225
generic mesh with rerefence counting
Definition geo.h:1089
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
std::ostream & os()
Definition diststream.h:247
the finite element space
Definition space.h:382
#define error_macro(message)
Definition dis_macros.h:49
#define warning_macro(message)
Definition dis_macros.h:53
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 vtk
size_type tensor_index(valued_type valued_tag, coordinate_type sys_coord, size_type i, size_type j)
This file is part of Rheolef.
odiststream & geo_put_vtk(odiststream &ops, const geo_basic< T, sequential > &omega, const basis_basic< T > &my_numb, const disarray< point_basic< T >, sequential > &my_node, bool append_data, size_t subgeo_dim)
odiststream & put_vtk_tensor_values(odiststream &ods, const field_basic< T, sequential > &tau_h, std::string name, bool put_header)
odiststream & put_vtk_scalar_values(odiststream &ods, const field_basic< T, sequential > &uh, std::string name, bool put_header)
template odiststream & field_put_vtk< Float >(odiststream &, const field_basic< Float, sequential > &, std::string, bool)
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
void put_header(odiststream &out, const branch_basic< T, sequential > &b)
Definition branch.cc:141
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
odiststream & put_vtk_vector_values(odiststream &ods, const field_basic< T, sequential > &uh, std::string name, bool put_header)
STL namespace.