ElementFileWriter.cc 11.1 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
namespace AMDiS {

  ElementFileWriter::ElementFileWriter(const std::string& name_, 
				       const FiniteElemSpace *feSpace_,
Thomas Witkowski's avatar
Thomas Witkowski committed
11
				       std::map<int, double> &mapvec)
12
13
14
15
16
17
18
19
20
21
22
23
    : 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),
Thomas Witkowski's avatar
Thomas Witkowski committed
24
      mesh(feSpace_->getMesh()),
25
      feSpace(feSpace_),
Thomas Witkowski's avatar
Thomas Witkowski committed
26
      vec(mapvec)
27
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
28
29
30
31
32
33
34
35
36
37
38
39
    if (name != "") {
      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);
      GET_PARAMETER(0, name + "->output->write every i-th timestep", "%d", &tsModulo);
    }
40
  }
41

42
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
45
46
  {
    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
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
98
99
100
101
102
103
104
  void ElementFileWriter::writeFile(std::map<int, double> &vec,
				    const FiniteElemSpace *feSpace,
				    const std::string& filename)
  {
    ElementFileWriter efw("", feSpace, vec);
    efw.writeVtkValues(filename);
  }
  

Thomas Witkowski's avatar
Thomas Witkowski committed
105
  void ElementFileWriter::writeTecPlotValues(const std::string &filename)
106
107
  {
    FUNCNAME("ElementFileWriter::writeTecPlotValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    std::ofstream fout(filename.c_str());
109
110
111

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

113
114
    int dim = mesh->getDim();
    double val;
115
116


117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
    // === 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;
145

146
147
148
149
150
151
152
153
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);

    while (elInfo) {
      // Get element value.
      val = vec[elInfo->getElement()->getIndex()];
154
    
155
156
      // Write coordinates of all element vertices and element value.
      for (int i = 0; i <= dim; ++i) {
157
      
158
159
160
161
	for (int j = 0; j < dim; ++j) {
	  fout << elInfo->getCoord(i)[j] << " ";
	}
	fout << val << "\n";
162
163
      }
    
164
165
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
166
167

  
168
169
170
171
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
    fout << "\n";
172

173
174
175
176
177
178
    for (int i = 0; i<numLeaves; ++i) {
      for (int j = 0; j<=dim; ++j) {
	++vertCntr;
	fout << vertCntr << " ";
      }
      fout << "\n";
179
    }
180
181
182


    fout.close();
183
184
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
185
  void ElementFileWriter::writeMeshDatValues(const std::string &filename, double time)
186
187
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
188
    std::ofstream fout(filename.c_str());
189

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

192
193
    int dim = mesh->getDim();
    double val;
194

195
196
197
198
199
200
201
    // === 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";
202
203


204
205
206
    // === Write vertex coordinates (every vertex for every element). ===
    fout << "vertex coordinates:\n";
    TraverseStack stack;
207

208
209
210
211
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
212

213
    while (elInfo) {
214

215
216
217
218
219
220
221
222
223
224
      // 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);
    }
225
226


227
228
229
230
231
232
233
234
235
    // === 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;
236
237
238
239
240
      }
      fout << "\n";
    }


241
    // === Write values. ===
242

243
    // Write values header.
244
    fout << "\n";
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
    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
273
274


275
276
277
278
    // Write values trailor.
    fout << "\n";
    fout << "interpolation values: " << name.c_str() << "\n\n\n";
    fout << "element interpolation points: " << name.c_str() << "\n";
279
280


281
282
    fout.close();
  }
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
  void ElementFileWriter::writeVtkValues(const std::string &filename)
285
286
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
287
    std::ofstream fout(filename.c_str());
288
289
290
291
292
293
294
295

    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;

296
297
298
    fout << "<?xml version=\"1.0\"?>" << std::endl;
    fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
    fout << "  <UnstructuredGrid>" << std::endl;
299
    fout << "    <Piece NumberOfPoints=\"" << (dim + 1) * mesh->getNumberOfLeaves() 
300
301
302
	 << "\" NumberOfCells=\"" <<  nElements << "\">" << std::endl;
    fout << "      <Points>" << std::endl;
    fout << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326


    // === 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);
    }
327

328
329
330
    fout << "        </DataArray>" << std::endl;
    fout << "      </Points>" << std::endl;
    fout << "      <Cells>" << std::endl;
331

332
    fout << "        <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
333
    for (int i = 0; i < nElements; i++) {
334
      fout << " " << (i + 1) * vertices << std::endl;
335
    }
336
    fout << "        </DataArray>" << std::endl;
337
338


339
    fout << "        <DataArray type=\"UInt8\" Name=\"types\">" << std::endl;
340
341
342
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
343
	fout << " 3" << std::endl;
344
345
	break;
      case 3: 
346
	fout << " 5" << std::endl;
347
348
	break;
      case 4:
349
	fout << " 10" << std::endl;
350
351
352
353
354
	break;		    
      default:
	break;
      }	   
    }
355
    fout << "        </DataArray>" << std::endl;
356

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

    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()];
382
   
383
384
385
386
387
388
389
      // 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);
    } 
390

391
    fout << "        </DataArray>" << std::endl;
392
    
393

394
395
396
397
    fout << "      </PointData>" << std::endl;
    fout << "    </Piece>" << std::endl;
    fout << "  </UnstructuredGrid>" << std::endl;
    fout << "</VTKFile>" << std::endl;
398

399
400
    fout.close();
  }
401
}