ElementFileWriter.cc 10.8 KB
Newer Older
1
2
3
4
5
6
#include "ElementFileWriter.h"
#include "BasisFunction.h"
#include "Parameters.h"
#include "Traverse.h"
#include "AdaptInfo.h"

7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
namespace AMDiS {

  ElementFileWriter::ElementFileWriter(const std::string& name_, 
				       Mesh *mesh_,
				       const FiniteElemSpace *feSpace_,
				       std::map<int, double> &vec_)
    : name(name_),
      tecplotExt(".plt"),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
      writeTecPlotFormat(0),
      writeAMDiSFormat(0),
      writeVtkFormat(0),
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
      mesh(mesh_),
      feSpace(feSpace_),
      vec(vec_)
  {
    GET_PARAMETER(0, name + "->output->filename", &filename);
    GET_PARAMETER(0, name + "->output->TecPlot format", "%d", 
		  &writeTecPlotFormat);
    GET_PARAMETER(0, name + "->output->TecPlot ext", &tecplotExt);
    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);
39
40
    GET_PARAMETER(0, name + "->output->write every i-th timestep", "%d", 
		  &tsModulo);
41
  }
42

43
44
45
46
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
				     int level, Flag traverseFlag, bool (*writeElem)(ElInfo*))
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
47

48
49
    timestepNumber++;
    timestepNumber %= tsModulo;
50

51
52
    if ((timestepNumber != 0) && !force) 
      return;
53

54
    std::string fn = filename;
55

56
57
58
59
    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");
60
    
61
      char formatStr[9];
62

63
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
64

65
      char timeStr[20];
66

67
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
68

69
70
      fn += timeStr;
    }
71
72


73
74
75
76
77
78
    if (writeTecPlotFormat) {
      TEST_EXIT(mesh)("no mesh\n");

      writeTecPlotValues(const_cast<char*>((fn + tecplotExt).c_str()));
      MSG("TecPlot file written to %s\n", (fn + tecplotExt).c_str());
    }
79

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

88
89
    if (writeVtkFormat) {
      TEST_EXIT(mesh)("no mesh\n");
90

91
92
93
      writeVtkValues(const_cast<char*>( (fn + vtkExt).c_str()));
      MSG("VTK file written to %s\n", (fn + vtkExt).c_str());		     
    }
94
95
  }

96
97
98
99
100
101
102
  void ElementFileWriter::writeTecPlotValues(const char* filename)
  {
    FUNCNAME("ElementFileWriter::writeTecPlotValues()");
    ::std::ofstream fout(filename);

    TEST_EXIT(fout)("Could not open file %s !\n", filename);
    fout.setf(std::ios::scientific,std::ios::floatfield);
103

104
105
    int dim = mesh->getDim();
    double val;
106
107


108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
    // === Write header. ===
    fout << "TITLE = \"" << name.c_str() << "\"\n";
    fout << "VARIABLES = ";
    switch (dim) {
    case 2: fout << "\"x\",\"y\"";
      break;
    case 3: fout << "\"x\",\"y\",\"z\"";
      break;
    default: ERROR_EXIT("illegal dimension !\n");
      break;
    }
    fout << ",\"" << name.c_str() << "\"\n";
    fout << "ZONE T=\"" << name.c_str() << "\""
	 << ", N=" << 3*mesh->getNumberOfLeaves() 
	 << ", E=" << mesh->getNumberOfLeaves() 
	 << ", F=FEPOINT, ";
    switch (dim) {
    case 2: fout << "ET=TRIANGLE\n\n";
      break;
    case 3: fout << "ET=TETRAHEDRON\n\n";
      break;
    default: ERROR_EXIT("illegal dimension !\n");
      break;
    }


    // === Write vertex coordinates and values (for each element !). ===
    TraverseStack stack;
136

137
138
139
140
141
142
143
144
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);

    while (elInfo) {
      // Get element value.
      val = vec[elInfo->getElement()->getIndex()];
145
    
146
147
      // Write coordinates of all element vertices and element value.
      for (int i = 0; i <= dim; ++i) {
148
      
149
150
151
152
	for (int j = 0; j < dim; ++j) {
	  fout << elInfo->getCoord(i)[j] << " ";
	}
	fout << val << "\n";
153
154
      }
    
155
156
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
157
158

  
159
160
161
162
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
    fout << "\n";
163

164
165
166
167
168
169
    for (int i = 0; i<numLeaves; ++i) {
      for (int j = 0; j<=dim; ++j) {
	++vertCntr;
	fout << vertCntr << " ";
      }
      fout << "\n";
170
    }
171
172
173


    fout.close();
174
175
  }

176
177
178
179
  void ElementFileWriter::writeMeshDatValues(const char* filename, double time)
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
    ::std::ofstream fout(filename);
180

181
    TEST_EXIT(fout)("Could not open file %s !\n", filename);
182

183
184
    int dim = mesh->getDim();
    double val;
185

186
187
188
189
190
191
192
    // === Write header. ===
    fout << "mesh name: " << mesh->getName().c_str() << "\n\n";
    fout << "time: " << time << "\n\n"; 
    fout << "DIM: " << dim << "\n";
    fout << "DIM_OF_WORLD: " << Global::getGeo(WORLD) << "\n\n";
    fout << "number of vertices: " << (dim+1)*mesh->getNumberOfLeaves() << "\n";
    fout << "number of elements: " << mesh->getNumberOfLeaves() << "\n\n";
193
194


195
196
197
    // === Write vertex coordinates (every vertex for every element). ===
    fout << "vertex coordinates:\n";
    TraverseStack stack;
198

199
200
201
202
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
203

204
    while (elInfo) {
205

206
207
208
209
210
211
212
213
214
215
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; ++i) {     
	for (int j = 0; j < dim; ++j) {
	  fout << elInfo->getCoord(i)[j] << " ";
	}
	fout << "\n";
      }
    
      elInfo = stack.traverseNext(elInfo);
    }
216
217


218
219
220
221
222
223
224
225
226
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
    fout << "\n";
    fout << "element vertices:\n";
    for (int i = 0; i < numLeaves; ++i) {
      for (int j = 0; j <= dim; ++j) {
	fout << vertCntr << " ";
	++vertCntr;
227
228
229
230
231
      }
      fout << "\n";
    }


232
    // === Write values. ===
233

234
    // Write values header.
235
    fout << "\n";
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
    fout << "number of values: 1\n\n";
    fout << "value description: " << name.c_str() << "\n";
    fout << "number of interpolation points: 0" << "\n";
    fout << "type: scalar" << "\n";
    fout << "interpolation type: lagrange" << "\n";
    fout << "interpolation degree: 1" << "\n";
    fout << "end of description: " << name.c_str() << "\n\n";

    // Write values.
    fout << "vertex values: " << name.c_str() << "\n";

    fout.setf(std::ios::scientific,std::ios::floatfield);

    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) {      
	fout << val << "\n";
      }
    
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
264
265


266
267
268
269
    // Write values trailor.
    fout << "\n";
    fout << "interpolation values: " << name.c_str() << "\n\n\n";
    fout << "element interpolation points: " << name.c_str() << "\n";
270
271


272
273
    fout.close();
  }
274

275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
  void ElementFileWriter::writeVtkValues(const char* filename)
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");
    ::std::ofstream fout(filename);

    TEST_EXIT(fout)("Could not open file %s !\n", filename);

    int dim = mesh->getDim();
    int vertices = mesh->getGeo(VERTEX);
    int nElements = mesh->getNumberOfLeaves();
    double val;

    fout << "<?xml version=\"1.0\"?>" << ::std::endl;
    fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << ::std::endl;
    fout << "  <UnstructuredGrid>" << ::std::endl;
    fout << "    <Piece NumberOfPoints=\"" << (dim + 1) * mesh->getNumberOfLeaves() 
	 << "\" NumberOfCells=\"" <<  nElements << "\">" << ::std::endl;
    fout << "      <Points>" << ::std::endl;
    fout << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << ::std::endl;


    // === Write vertex coordinates (every vertex for every element). ===
    TraverseStack stack;

    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);

    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {     
	for (int j = 0; j < dim; j++) {
	  fout << elInfo->getCoord(i)[j] << " ";
	}
	for (int j = dim; j < 3; j++) {
	  fout << "0.0";
	}
	fout << "\n";
      }
   
      elInfo = stack.traverseNext(elInfo);
    }
318

319
320
321
    fout << "        </DataArray>" << ::std::endl;
    fout << "      </Points>" << ::std::endl;
    fout << "      <Cells>" << ::std::endl;
322

323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
    fout << "        <DataArray type=\"Int32\" Name=\"offsets\">" << ::std::endl;
    for (int i = 0; i < nElements; i++) {
      fout << " " << (i + 1) * vertices << ::std::endl;
    }
    fout << "        </DataArray>" << ::std::endl;


    fout << "        <DataArray type=\"UInt8\" Name=\"types\">" << ::std::endl;
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
	fout << " 3" << ::std::endl;
	break;
      case 3: 
	fout << " 5" << ::std::endl;
	break;
      case 4:
	fout << " 10" << ::std::endl;
	break;		    
      default:
	break;
      }	   
    }
    fout << "        </DataArray>" << ::std::endl;

    fout << "        <DataArray type=\"Int32\" Name=\"connectivity\">" << ::std::endl;
    int vertCntr = 0;
    for (int i = 0; i < nElements; ++i) {
      for (int j = 0; j <= dim; ++j) {
	fout << vertCntr << " ";
	++vertCntr;
      }
      fout << ::std::endl;
    }
    fout << "        </DataArray>" << ::std::endl;
358
    
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    fout << "      </Cells>" << ::std::endl;
    fout << "      <PointData>" << ::std::endl;
    fout << "        <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">" << ::std::endl;

    fout.setf(std::ios::scientific,std::ios::floatfield);

    elInfo = stack.traverseFirst(mesh,
				 -1, 
				 Mesh::CALL_LEAF_EL |
				 Mesh::FILL_COORDS);
    int vc = 0;
    while (elInfo) {
      // Get element value.
      val = vec[elInfo->getElement()->getIndex()];
373
   
374
375
376
377
378
379
380
      // Write value for each vertex of each element.
      for (int i = 0; i <= dim; i++) {      
	fout << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
      }
      vc++;
      elInfo = stack.traverseNext(elInfo);
    } 
381

382
383
    fout << "        </DataArray>" << ::std::endl;
    
384

385
386
387
388
    fout << "      </PointData>" << ::std::endl;
    fout << "    </Piece>" << ::std::endl;
    fout << "  </UnstructuredGrid>" << ::std::endl;
    fout << "</VTKFile>" << ::std::endl;
389

390
391
    fout.close();
  }
392
}