ElementFileWriter.cc 12 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
39
40
41
42
43
44
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);
    GET_PARAMETER(0, name + "->output->write every i-th timestep", "%d", 
		  &tsModulo);

    //   TEST_EXIT(vec.size() == mesh->getNumberOfLeaves())
    //     ("illegal size !\n");
  }
45

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

51
52
    timestepNumber++;
    timestepNumber %= tsModulo;
53

54
55
    if ((timestepNumber != 0) && !force) 
      return;
56

57
    std::string fn = filename;
58

59
60
61
62
    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");
63
    
64
      char formatStr[9];
65

66
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
67

68
      char timeStr[20];
69

70
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
71

72
73
      fn += timeStr;
    }
74
75


76
77
78
79
80
81
    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());
    }
82

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

91
92
    if (writeVtkFormat) {
      TEST_EXIT(mesh)("no mesh\n");
93

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

99
100
101
102
103
104
105
  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);
106

107
108
    int dim = mesh->getDim();
    double val;
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
136
137
138
    // === 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;
139

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

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

  
162
163
164
165
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
    fout << "\n";
166

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


    fout.close();
177
178
  }

179
180
181
182
  void ElementFileWriter::writeMeshDatValues(const char* filename, double time)
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
    ::std::ofstream fout(filename);
183

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

186
187
    int dim = mesh->getDim();
    double val;
188

189
190
191
192
193
194
195
    // === 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";
196
197


198
199
200
    // === Write vertex coordinates (every vertex for every element). ===
    fout << "vertex coordinates:\n";
    TraverseStack stack;
201

202
203
204
205
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
206

207
    while (elInfo) {
208

209
210
211
212
213
214
215
216
217
218
      // 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);
    }
219
220


221
222
223
224
225
226
227
228
229
    // === 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;
230
231
232
233
234
      }
      fout << "\n";
    }


235
    // === Write values. ===
236

237
    // Write values header.
238
    fout << "\n";
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
264
265
266
    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
267
268


269
270
271
272
    // Write values trailor.
    fout << "\n";
    fout << "interpolation values: " << name.c_str() << "\n\n\n";
    fout << "element interpolation points: " << name.c_str() << "\n";
273
274


275
276
    fout.close();
  }
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
318
319
320
  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);
    }
321

322
323
324
    fout << "        </DataArray>" << ::std::endl;
    fout << "      </Points>" << ::std::endl;
    fout << "      <Cells>" << ::std::endl;
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
358
359
360
    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;
361
    
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
    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()];

      /*      if ((elInfo->getCoord(0)[0] == 80.0) &&
	  (elInfo->getCoord(0)[1] == 1.25) &&
	  (elInfo->getCoord(1)[0] == 79.9609375) &&
	  (elInfo->getCoord(1)[1] == 1.2109375) &&
	  (elInfo->getCoord(2)[0] == 80.0) &&
	  (elInfo->getCoord(2)[1] == 1.2109375))
	{
	  ::std::cout << "------" << ::std::endl;
	  ::std::cout << "val = " << val << "  I = " << elInfo->getElement()->getIndex() << ::std::endl;
	  ::std::cout << "PARENT-I = " << elInfo->getParent()->getIndex() << ::std::endl;
	  ::std::cout << "vc = " << vc << ::std::endl;
	  ::std::cout.precision(10);
	  for (int i = 0; i < 3; i++) {
	    ::std::cout << elInfo->getCoord(i)[0] << "/" << elInfo->getCoord(i)[1] << ::std::endl;
	  }	  	  
	}
      */
      
      if (val > 0.3) {
	::std::cout.precision(10);
	::std::cout << "------" << ::std::endl;
	::std::cout << "val = " << val << "  I = " << elInfo->getElement()->getIndex() << ::std::endl;
	::std::cout << "vc = " << vc << ::std::endl;
	for (int i = 0; i < 3; i++) {
	  ::std::cout << elInfo->getCoord(i)[0] << "/" << elInfo->getCoord(i)[1] << ::std::endl;
	}	  	  
      }
404
405
      
    
406
407
408
409
410
411
412
      // 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);
    } 
413

414
415
    fout << "        </DataArray>" << ::std::endl;
    
416

417
418
419
420
    fout << "      </PointData>" << ::std::endl;
    fout << "    </Piece>" << ::std::endl;
    fout << "  </UnstructuredGrid>" << ::std::endl;
    fout << "</VTKFile>" << ::std::endl;
421

422
423
    fout.close();
  }
424
}