ElementFileWriter.cc 9.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/convenience.hpp>

18
19
20
21
22
23
#include "ElementFileWriter.h"
#include "BasisFunction.h"
#include "Parameters.h"
#include "Traverse.h"
#include "AdaptInfo.h"

24
25
namespace AMDiS {

26
27
28
  using namespace std;

  ElementFileWriter::ElementFileWriter(string name_, 
29
				       Mesh *mesh_,
30
				       map<int, double> &mapvec)
31
32
33
34
35
36
37
38
39
40
    : name(name_),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
      writeAMDiSFormat(0),
      writeVtkFormat(0),
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
41
      mesh(mesh_),
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      vec(mapvec)
43
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
44
45
46
47
48
49
50
51
52
53
    if (name != "") {
      GET_PARAMETER(0, name + "->output->filename", &filename);
      GET_PARAMETER(0, name + "->output->AMDiS format", "%d", &writeAMDiSFormat);
      GET_PARAMETER(0, name + "->output->AMDiS mesh-dat ext", &amdisMeshDatExt);
      GET_PARAMETER(0, name + "->output->ParaView format", "%d", &writeVtkFormat);
      GET_PARAMETER(0, name + "->output->append index", "%d", &appendIndex);
      GET_PARAMETER(0, name + "->output->index length", "%d", &indexLength);
      GET_PARAMETER(0, name + "->output->index decimals", "%d", &indexDecimals);
      GET_PARAMETER(0, name + "->output->write every i-th timestep", "%d", &tsModulo);
    }
54
  }
55

56

57
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
60
61
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
62

63
64
    timestepNumber++;
    timestepNumber %= tsModulo;
65

66
    if (timestepNumber != 0 && !force) 
67
      return;
68

69
    string fn = filename;
70

71
72
73
74
    if (appendIndex) {
      TEST_EXIT(indexLength <= 99)("index lenght > 99\n");
      TEST_EXIT(indexDecimals <= 97)("index decimals > 97\n");
      TEST_EXIT(indexDecimals < indexLength)("index length <= index decimals\n");
75
    
76
77
      char formatStr[9];
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
78

79
80
      char timeStr[20];
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
81

82
83
      fn += timeStr;
    }
84

85
86
    if (writeAMDiSFormat) {
      TEST_EXIT(mesh)("no mesh\n");
87
    
88
      writeMeshDatValues(const_cast<char*>((fn + amdisMeshDatExt).c_str()),
89
90
91
			 adaptInfo ? adaptInfo->getTime() : 0.0);
      MSG("MeshDat file written to %s\n", (fn + amdisMeshDatExt).c_str());
    }
92

93
94
    if (writeVtkFormat) {
      TEST_EXIT(mesh)("no mesh\n");
95

96
      writeVtkValues(const_cast<char*>((fn + vtkExt).c_str()));
97
98
      MSG("VTK file written to %s\n", (fn + vtkExt).c_str());		     
    }
99
100
  }

101

102
  void ElementFileWriter::writeFile(map<int, double> &vec,
103
				    Mesh *mesh,
104
				    string filename,
105
				    int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
106
  {
107
    ElementFileWriter efw("", mesh, vec);
108
    efw.writeVtkValues(filename, level);
Thomas Witkowski's avatar
Thomas Witkowski committed
109
110
111
  }
  

112
  void ElementFileWriter::writeMeshDatValues(string filename, double time)
113
114
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
115

116
117
118
119
120
121
122
123
124
    boost::iostreams::filtering_ostream file;
    {
      //boost::iostreams seems not to truncate the file
      ofstream swapfile(filename.c_str(), ios::out | ios::trunc);
      TEST_EXIT(swapfile.is_open())
	("Cannot open file %s for writing!\n", name.c_str());
      swapfile.close();
    }
    file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc));
125

126
127
    int dim = mesh->getDim();
    double val;
128

129
    // === Write header. ===
130
131
132
133
134
135
    file << "mesh name: " << mesh->getName().c_str() << "\n\n";
    file << "time: " << time << "\n\n"; 
    file << "DIM: " << dim << "\n";
    file << "DIM_OF_WORLD: " << Global::getGeo(WORLD) << "\n\n";
    file << "number of vertices: " << (dim+1)*mesh->getNumberOfLeaves() << "\n";
    file << "number of elements: " << mesh->getNumberOfLeaves() << "\n\n";
136
137


138
    // === Write vertex coordinates (every vertex for every element). ===
139
    file << "vertex coordinates:\n";
140
    TraverseStack stack;
141

142
143
144
145
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
146

147
    while (elInfo) {
148

149
150
151
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; ++i) {     
	for (int j = 0; j < dim; ++j) {
152
	  file << elInfo->getCoord(i)[j] << " ";
153
	}
154
	file << "\n";
155
156
157
158
      }
    
      elInfo = stack.traverseNext(elInfo);
    }
159
160


161
162
163
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
164
165
    file << "\n";
    file << "element vertices:\n";
166
167
    for (int i = 0; i < numLeaves; ++i) {
      for (int j = 0; j <= dim; ++j) {
168
	file << vertCntr << " ";
169
	++vertCntr;
170
      }
171
      file << "\n";
172
173
174
    }


175
    // === Write values. ===
176

177
    // Write values header.
178
179
180
181
182
183
184
185
    file << "\n";
    file << "number of values: 1\n\n";
    file << "value description: " << name.c_str() << "\n";
    file << "number of interpolation points: 0" << "\n";
    file << "type: scalar" << "\n";
    file << "interpolation type: lagrange" << "\n";
    file << "interpolation degree: 1" << "\n";
    file << "end of description: " << name.c_str() << "\n\n";
186
187

    // Write values.
188
    file << "vertex values: " << name.c_str() << "\n";
189

190
    file.setf(ios::scientific,ios::floatfield);
191
192
193
194
195
196
197
198
199
200
201

    elInfo = stack.traverseFirst(mesh,
				 -1, 
				 Mesh::CALL_LEAF_EL);

    while (elInfo) {
      // Get element value.
      val = vec[elInfo->getElement()->getIndex()];
    
      // Write value for each vertex of each element.
      for (int i = 0; i <= dim; ++i) {      
202
	file << val << "\n";
203
204
205
206
      }
    
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
207
208


209
    // Write values trailor.
210
211
212
    file << "\n";
    file << "interpolation values: " << name.c_str() << "\n\n\n";
    file << "element interpolation points: " << name.c_str() << "\n";
213
  }
214

215

216
  void ElementFileWriter::writeVtkValues(string filename, int level)
217
218
219
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");

220
221
222
223
224
225
226
227
228
    boost::iostreams::filtering_ostream file;
    {
      //boost::iostreams seems not to truncate the file
      ofstream swapfile(filename.c_str(), ios::out | ios::trunc);
      TEST_EXIT(swapfile.is_open())
	("Cannot open file %s for writing!\n", name.c_str());
      swapfile.close();
    }
    file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc));
229
230

    int dim = mesh->getDim();
231
    int dimOfWorld = Global::getGeo(WORLD);
232
233
    int vertices = mesh->getGeo(VERTEX);
    int nElements = mesh->getNumberOfLeaves();
234
235
236
237
238
239
240
241
242
243
244
245
    Flag traverseFlag = level == -1 ? Mesh::CALL_LEAF_EL : Mesh::CALL_EL_LEVEL;

    if (level != -1) {
      nElements = 0;

      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, level, Mesh::CALL_EL_LEVEL);      
      while (elInfo) {
	nElements++;
	elInfo = stack.traverseNext(elInfo);
      }
    }
246

247
248
249
250
251
252
253
    file << "<?xml version=\"1.0\"?>\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "  <UnstructuredGrid>\n";
    file << "    <Piece NumberOfPoints=\"" << (dim + 1) * nElements
	 << "\" NumberOfCells=\"" <<  nElements << "\">\n";
    file << "      <Points>\n";
    file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
254
255
256
257


    // === Write vertex coordinates (every vertex for every element). ===
    TraverseStack stack;
258
    ElInfo *elInfo = 
259
      stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
260
261
262
263

    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {     
264
	for (int j = 0; j < dimOfWorld; j++) 
265
	  file << elInfo->getCoord(i)[j] << " ";
266
267
	
	for (int j = dimOfWorld; j < 3; j++) 
268
	  file << "0.0";
269
	
270
	file << "\n";
271
272
273
274
      }
   
      elInfo = stack.traverseNext(elInfo);
    }
275

276
277
278
    file << "        </DataArray>\n";
    file << "      </Points>\n";
    file << "      <Cells>\n";
279

280
    file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
281
    for (int i = 0; i < nElements; i++) 
282
283
      file << " " << (i + 1) * vertices << "\n";    
    file << "        </DataArray>\n";
284
285


286
    file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";
287
288
289
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
290
	file << " 3\n";
291
292
	break;
      case 3: 
293
	file << " 5\n";
294
295
	break;
      case 4:
296
	file << " 10\n";
297
298
299
300
301
	break;		    
      default:
	break;
      }	   
    }
302
    file << "        </DataArray>\n";
303

304
    file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
305
306
307
    int vertCntr = 0;
    for (int i = 0; i < nElements; ++i) {
      for (int j = 0; j <= dim; ++j) {
308
	file << vertCntr << " ";
309
310
	++vertCntr;
      }
311
      file << "\n";
312
    }
313
    file << "        </DataArray>\n";
314
    
315
316
317
    file << "      </Cells>\n";
    file << "      <PointData>\n";
    file << "        <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">\n";
318

319
    file.setf(ios::scientific,ios::floatfield);
320

321
    elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
322
323
324
    int vc = 0;
    while (elInfo) {
      // Get element value.
325
      double val = vec[elInfo->getElement()->getIndex()];
326
   
327
      // Write value for each vertex of each element.
328
      for (int i = 0; i <= dim; i++) 
329
	file << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
330
      
331
332
333
      vc++;
      elInfo = stack.traverseNext(elInfo);
    } 
334

335
    file << "        </DataArray>\n";
336
    
337

338
339
340
341
    file << "      </PointData>\n";
    file << "    </Piece>\n";
    file << "  </UnstructuredGrid>\n";
    file << "</VTKFile>\n";
342
  }
343
}