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

9
  ElementFileWriter::ElementFileWriter(std::string name_, 
10
				       Mesh *mesh_,
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),
24
      mesh(mesh_),
Thomas Witkowski's avatar
Thomas Witkowski committed
25
      vec(mapvec)
26
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
27 28 29 30 31 32 33 34 35 36 37 38
    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);
    }
39
  }
40

41
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
42 43
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
44 45
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
46

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

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

53
    std::string fn = filename;
54

55 56 57 58
    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");
59
    
60 61
      char formatStr[9];
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
62

63 64
      char timeStr[20];
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
65

66 67
      fn += timeStr;
    }
68 69


70 71 72 73 74 75
    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());
    }
76

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

85 86
    if (writeVtkFormat) {
      TEST_EXIT(mesh)("no mesh\n");
87

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

93

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

103
  void ElementFileWriter::writeTecPlotValues(std::string filename)
104 105
  {
    FUNCNAME("ElementFileWriter::writeTecPlotValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
106
    std::ofstream fout(filename.c_str());
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108
    TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
109
    fout.setf(std::ios::scientific,std::ios::floatfield);
110

111 112
    int dim = mesh->getDim();
    double val;
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 139 140 141 142
    // === 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;
143

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

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

  
165 166 167 168
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
    fout << "\n";
169

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


    fout.close();
180 181
  }

182

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

Thomas Witkowski's avatar
Thomas Witkowski committed
188
    TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
189

190 191
    int dim = mesh->getDim();
    double val;
192

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


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

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

211
    while (elInfo) {
212

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


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


239
    // === Write values. ===
240

241
    // Write values header.
242
    fout << "\n";
243 244 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
    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
271 272


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


279 280
    fout.close();
  }
281

282

283
  void ElementFileWriter::writeVtkValues(std::string filename)
284 285
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    std::ofstream fout(filename.c_str());
287

Thomas Witkowski's avatar
Thomas Witkowski committed
288
    TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
289 290

    int dim = mesh->getDim();
291
    int dimOfWorld = Global::getGeo(WORLD);
292 293 294
    int vertices = mesh->getGeo(VERTEX);
    int nElements = mesh->getNumberOfLeaves();

295 296 297
    fout << "<?xml version=\"1.0\"?>" << std::endl;
    fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
    fout << "  <UnstructuredGrid>" << std::endl;
298
    fout << "    <Piece NumberOfPoints=\"" << (dim + 1) * nElements
299 300 301
	 << "\" NumberOfCells=\"" <<  nElements << "\">" << std::endl;
    fout << "      <Points>" << std::endl;
    fout << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
302 303 304 305


    // === Write vertex coordinates (every vertex for every element). ===
    TraverseStack stack;
306 307
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
308 309 310 311

    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {     
312
	for (int j = 0; j < dimOfWorld; j++) 
313
	  fout << elInfo->getCoord(i)[j] << " ";
314 315
	
	for (int j = dimOfWorld; j < 3; j++) 
316
	  fout << "0.0";
317
	
318 319 320 321 322
	fout << "\n";
      }
   
      elInfo = stack.traverseNext(elInfo);
    }
323

324 325 326
    fout << "        </DataArray>" << std::endl;
    fout << "      </Points>" << std::endl;
    fout << "      <Cells>" << std::endl;
327

328
    fout << "        <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
329 330
    for (int i = 0; i < nElements; i++) 
      fout << " " << (i + 1) * vertices << std::endl;    
331
    fout << "        </DataArray>" << std::endl;
332 333


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

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

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

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

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

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

391 392
    fout.close();
  }
393
}