1 /*
2 * Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 using namespace std;
19
20 #include "vtk_file_writer.h"
21
22 #include <vtkRectilinearGrid.h>
23 #include <vtkRectilinearGridWriter.h>
24 #include <vtkXMLRectilinearGridWriter.h>
25 #include <vtkStructuredGrid.h>
26 #include <vtkStructuredGridWriter.h>
27 #include <vtkXMLStructuredGridWriter.h>
28 #include <vtkZLibDataCompressor.h>
29 #include <vtkFloatArray.h>
30 #include <vtkDoubleArray.h>
31 #include <vtkFieldData.h>
32 #include <vtkPointData.h>
33
34 #include <sstream>
35 #include <iomanip>
36
37
VTK_File_Writer(string filename,int meshType)38 VTK_File_Writer::VTK_File_Writer(string filename, int meshType)
39 {
40 SetFilename(filename);
41 m_MeshType = meshType;
42 m_NativeDump = false;
43 m_Binary = true;
44 m_Compress = true;
45 m_AppendMode = false;
46 m_ActiveTS = false;
47
48 if (m_MeshType==0) //cartesian mesh
49 m_GridData = vtkRectilinearGrid::New();
50 else if (m_MeshType==1) //cylindrical mesh
51 m_GridData = vtkStructuredGrid::New();
52 else
53 {
54 cerr << "VTK_File_Writer::VTK_File_Writer: Error, unknown mesh type: " << m_MeshType << endl;
55 m_GridData=NULL;
56 }
57 }
58
~VTK_File_Writer()59 VTK_File_Writer::~VTK_File_Writer()
60 {
61 if (m_GridData)
62 m_GridData->Delete();
63 m_GridData = NULL;
64 }
65
SetMeshLines(double const * const * lines,unsigned int const * count,double scaling)66 void VTK_File_Writer::SetMeshLines(double const* const* lines, unsigned int const* count, double scaling)
67 {
68 if (m_MeshType==0) //cartesian mesh
69 {
70 vtkRectilinearGrid* RectGrid = dynamic_cast<vtkRectilinearGrid*>(m_GridData);
71 if (RectGrid==NULL)
72 {
73 cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
74 exit(1);
75 }
76 RectGrid->SetDimensions(count[0],count[1],count[2]);
77 vtkDoubleArray *Coords[3];
78 for (int n=0;n<3;++n)
79 {
80 m_MeshLines[n].clear();
81 m_MeshLines[n].reserve(count[n]);
82 Coords[n] = vtkDoubleArray::New();
83 for (unsigned int i=0; i<count[n]; i++)
84 {
85 Coords[n]->InsertNextValue(lines[n][i]*scaling);
86 m_MeshLines[n].push_back(lines[n][i]*scaling);
87 }
88 }
89 RectGrid->SetXCoordinates(Coords[0]);
90 RectGrid->SetYCoordinates(Coords[1]);
91 RectGrid->SetZCoordinates(Coords[2]);
92 for (int n=0;n<3;++n)
93 Coords[n]->Delete();
94 }
95 else if (m_MeshType==1) //cylindrical mesh
96 {
97 vtkStructuredGrid* StructGrid = dynamic_cast<vtkStructuredGrid*>(m_GridData);
98 if (StructGrid==NULL)
99 {
100 cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
101 exit(1);
102 }
103
104 for (int n=0;n<3;++n)
105 {
106 m_MeshLines[n].clear();
107 m_MeshLines[n].reserve(count[n]);
108 double scale=1;
109 if (n!=1)
110 scale*=scaling;
111 for (unsigned int i=0; i<count[n]; i++)
112 m_MeshLines[n].push_back(lines[n][i]*scale);
113 }
114
115 StructGrid->SetDimensions(count[0],count[1],count[2]);
116 vtkPoints *points = vtkPoints::New();
117 points->SetNumberOfPoints(count[0]*count[1]*count[2]);
118 double r[3];
119 int id=0;
120 for (unsigned int k=0; k<count[2]; ++k)
121 for (unsigned int j=0; j<count[1]; ++j)
122 for (unsigned int i=0; i<count[0]; ++i)
123 {
124 r[0] = lines[0][i] * cos(lines[1][j]) * scaling;
125 r[1] = lines[0][i] * sin(lines[1][j]) * scaling;
126 r[2] = lines[2][k] * scaling;
127 points->SetPoint(id++,r);
128 }
129 StructGrid->SetPoints(points);
130 points->Delete();
131 }
132 else
133 {
134 cerr << "VTK_File_Writer::SetMeshLines: Error, unknown mesh type: " << m_MeshType << endl;
135 }
136 }
137
AddScalarField(string fieldname,double const * const * const * field)138 void VTK_File_Writer::AddScalarField(string fieldname, double const* const* const* field)
139 {
140 vtkDoubleArray* array = vtkDoubleArray::New();
141 array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
142 array->SetName(fieldname.c_str());
143 int id=0;
144 for (unsigned int k=0;k<m_MeshLines[2].size();++k)
145 {
146 for (unsigned int j=0;j<m_MeshLines[1].size();++j)
147 {
148 for (unsigned int i=0;i<m_MeshLines[0].size();++i)
149 {
150 array->SetTuple1(id++,field[i][j][k]);
151 }
152 }
153 }
154 m_GridData->GetPointData()->AddArray(array);
155 array->Delete();
156 }
157
AddScalarField(string fieldname,float const * const * const * field)158 void VTK_File_Writer::AddScalarField(string fieldname, float const* const* const* field)
159 {
160 vtkFloatArray* array = vtkFloatArray::New();
161 array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
162 array->SetName(fieldname.c_str());
163 int id=0;
164 for (unsigned int k=0;k<m_MeshLines[2].size();++k)
165 {
166 for (unsigned int j=0;j<m_MeshLines[1].size();++j)
167 {
168 for (unsigned int i=0;i<m_MeshLines[0].size();++i)
169 {
170 array->SetTuple1(id++,field[i][j][k]);
171 }
172 }
173 }
174 m_GridData->GetPointData()->AddArray(array);
175 array->Delete();
176 }
177
AddVectorField(string fieldname,double const * const * const * const * field)178 void VTK_File_Writer::AddVectorField(string fieldname, double const* const* const* const* field)
179 {
180 vtkDoubleArray* array = vtkDoubleArray::New();
181 array->SetNumberOfComponents(3);
182 array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
183 array->SetName(fieldname.c_str());
184 int id=0;
185 double out[3];
186 for (unsigned int k=0;k<m_MeshLines[2].size();++k)
187 {
188 for (unsigned int j=0;j<m_MeshLines[1].size();++j)
189 {
190 double cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
191 double sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
192 for (unsigned int i=0;i<m_MeshLines[0].size();++i)
193 {
194 if ((m_MeshType==0) || (m_NativeDump))
195 array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
196 else
197 {
198 out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
199 out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
200 out[2] = field[2][i][j][k];
201 array->SetTuple3(id++,out[0],out[1],out[2]);
202 }
203 }
204 }
205 }
206 m_GridData->GetPointData()->AddArray(array);
207 array->Delete();
208 }
209
AddVectorField(string fieldname,float const * const * const * const * field)210 void VTK_File_Writer::AddVectorField(string fieldname, float const* const* const* const* field)
211 {
212 vtkFloatArray* array = vtkFloatArray::New();
213 array->SetNumberOfComponents(3);
214 array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
215 array->SetName(fieldname.c_str());
216 int id=0;
217 float out[3];
218 for (unsigned int k=0;k<m_MeshLines[2].size();++k)
219 {
220 for (unsigned int j=0;j<m_MeshLines[1].size();++j)
221 {
222 float cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
223 float sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
224 for (unsigned int i=0;i<m_MeshLines[0].size();++i)
225 {
226 if ((m_MeshType==0) || (m_NativeDump))
227 array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
228 else
229 {
230 out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
231 out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
232 out[2] = field[2][i][j][k];
233 array->SetTuple3(id++,out[0],out[1],out[2]);
234 }
235 }
236 }
237 }
238 m_GridData->GetPointData()->AddArray(array);
239 array->Delete();
240 }
241
242
GetNumberOfFields() const243 int VTK_File_Writer::GetNumberOfFields() const
244 {
245 return m_GridData->GetPointData()->GetNumberOfArrays();
246 }
247
ClearAllFields()248 void VTK_File_Writer::ClearAllFields()
249 {
250 while (m_GridData->GetPointData()->GetNumberOfArrays()>0)
251 {
252 const char* name = m_GridData->GetPointData()->GetArrayName(0);
253 m_GridData->GetPointData()->RemoveArray(name);
254 }
255 }
256
Write()257 bool VTK_File_Writer::Write()
258 {
259 return WriteXML();
260 }
261
GetTimestepFilename(int pad_length) const262 string VTK_File_Writer::GetTimestepFilename(int pad_length) const
263 {
264 if (m_ActiveTS==false)
265 return m_filename;
266
267 stringstream ss;
268 ss << m_filename << "_" << std::setw( pad_length ) << std::setfill( '0' ) << m_timestep;
269
270 return ss.str();
271 }
272
273
WriteASCII()274 bool VTK_File_Writer::WriteASCII()
275 {
276 vtkDataWriter* writer = NULL;
277 if (m_MeshType==0) //cartesian mesh
278 writer = vtkRectilinearGridWriter::New();
279 else if (m_MeshType==1) //cylindrical mesh
280 writer = vtkStructuredGridWriter::New();
281 else
282 {
283 cerr << "VTK_File_Writer::WriteASCII: Error, unknown mesh type: " << m_MeshType << endl;
284 return false;
285 }
286
287 writer->SetHeader(m_header.c_str());
288 #if VTK_MAJOR_VERSION>=6
289 writer->SetInputData(m_GridData);
290 #else
291 writer->SetInput(m_GridData);
292 #endif
293
294 string filename = GetTimestepFilename() + ".vtk";
295 writer->SetFileName(filename.c_str());
296 if (m_Binary)
297 writer->SetFileTypeToBinary();
298 else
299 writer->SetFileTypeToASCII();
300
301 writer->Write();
302 writer->Delete();
303 return true;
304 }
305
WriteXML()306 bool VTK_File_Writer::WriteXML()
307 {
308 vtkXMLStructuredDataWriter* writer = NULL;
309 if (m_MeshType==0) //cartesian mesh
310 writer = vtkXMLRectilinearGridWriter::New();
311 else if (m_MeshType==1) //cylindrical mesh
312 writer = vtkXMLStructuredGridWriter::New();
313 else
314 {
315 cerr << "VTK_File_Writer::WriteXML: Error, unknown mesh type: " << m_MeshType << endl;
316 return false;
317 }
318
319 #if VTK_MAJOR_VERSION>=6
320 writer->SetInputData(m_GridData);
321 #else
322 writer->SetInput(m_GridData);
323 #endif
324
325 string filename = GetTimestepFilename() + "." + writer->GetDefaultFileExtension();
326 writer->SetFileName(filename.c_str());
327 if (m_Compress)
328 writer->SetCompressor(vtkZLibDataCompressor::New());
329 else
330 writer->SetCompressor(NULL);
331
332 if (m_Binary)
333 writer->SetDataModeToBinary();
334 else
335 writer->SetDataModeToAscii();
336
337 writer->Write();
338 writer->Delete();
339 return true;
340 }
341