ElementFileWriter.cc 13.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20 21


22 23 24 25
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/convenience.hpp>
26
#include <boost/lexical_cast.hpp>
27

28
#include "io/VtkWriter.h"
29 30
#include "ElementFileWriter.h"
#include "BasisFunction.h"
31
#include "Initfile.h"
32 33 34
#include "Traverse.h"
#include "AdaptInfo.h"

35
namespace AMDiS { namespace io {
36

37 38 39
  using namespace std;

  ElementFileWriter::ElementFileWriter(string name_, 
40
				       Mesh *mesh_,
41
				       map<int, double> &mapvec)
42 43 44
    : name(name_),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
45
      pvdExt(".pvd"),
46 47
      writeAMDiSFormat(0),
      writeVtkFormat(0),
48 49 50
      writeVtkVectorFormat(0),
      writeAs3dVector(true),
      writeParaViewAnimation(0),
51 52 53 54 55
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
56
      mesh(mesh_),
57
      vec(&mapvec),
58
      vecs(NULL)
59
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
60
    if (name != "") {
61 62 63 64
      Parameters::get(name + "->output->filename", filename);
      Parameters::get(name + "->output->AMDiS format", writeAMDiSFormat);
      Parameters::get(name + "->output->AMDiS mesh-dat ext", amdisMeshDatExt);
      Parameters::get(name + "->output->ParaView format", writeVtkFormat);
65
      Parameters::get(name + "->output->ParaView animation", writeParaViewAnimation);
66 67 68 69
      Parameters::get(name + "->output->append index", appendIndex);
      Parameters::get(name + "->output->index length", indexLength);
      Parameters::get(name + "->output->index decimals", indexDecimals);
      Parameters::get(name + "->output->write every i-th timestep", tsModulo);
Thomas Witkowski's avatar
Thomas Witkowski committed
70
    }
71
  }
72

73

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
  ElementFileWriter::ElementFileWriter(string name_,
				       Mesh *mesh_,
				       map<int, vector<double> > &mapvec)
    : name(name_),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
      pvdExt(".pvd"),
      writeAMDiSFormat(0),
      writeVtkFormat(0),
      writeVtkVectorFormat(0),
      writeAs3dVector(true),
      writeParaViewAnimation(0),
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
      mesh(mesh_),
92
      vec(NULL),
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
      vecs(&mapvec)
  {
    if (name != "") {
      Parameters::get(name + "->output->filename", filename);
      Parameters::get(name + "->output->AMDiS format", writeAMDiSFormat);
      Parameters::get(name + "->output->AMDiS mesh-dat ext", amdisMeshDatExt);
      Parameters::get(name + "->output->ParaView format", writeVtkFormat);
      Parameters::get(name + "->output->ParaView vector format", writeVtkVectorFormat);
      Parameters::get(name + "->output->write vector as 3d vector", writeAs3dVector);
      Parameters::get(name + "->output->ParaView animation", writeParaViewAnimation);
      Parameters::get(name + "->output->append index", appendIndex);
      Parameters::get(name + "->output->index length", indexLength);
      Parameters::get(name + "->output->index decimals", indexDecimals);
      Parameters::get(name + "->output->write every i-th timestep", tsModulo);
    }
  }
  

111
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
112 113
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
114 115
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
116

117 118
    timestepNumber++;
    timestepNumber %= tsModulo;
119

120
    if (timestepNumber != 0 && !force) 
121
      return;
122

123
    string fn = filename;
124

125 126 127 128
    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");
129
    
130 131
      char formatStr[9];
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
132

133 134
      char timeStr[20];
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
135

136 137
      fn += timeStr;
    }
138

139 140
    if (writeAMDiSFormat) {
      TEST_EXIT(mesh)("no mesh\n");
141
    
142
      writeMeshDatValues(const_cast<char*>((fn + amdisMeshDatExt).c_str()),
143 144 145
			 adaptInfo ? adaptInfo->getTime() : 0.0);
      MSG("MeshDat file written to %s\n", (fn + amdisMeshDatExt).c_str());
    }
146

147
    if (writeVtkFormat || writeVtkVectorFormat) {
148
      TEST_EXIT(mesh)("no mesh\n");
149

150
      writeVtkValues(fn, vtkExt, -1, writeVtkVectorFormat == 1);
151 152
      MSG("VTK file written to %s\n", (fn + vtkExt).c_str());		     
    }
153

154
    if (writeParaViewAnimation)
155
      VtkWriter::detail::updateAnimationFile(adaptInfo, 
156 157 158
				     (fn + vtkExt), 
				     &paraViewAnimationFrames, 
				     (filename + pvdExt));
159 160
  }

161

162
  void ElementFileWriter::writeFile(map<int, double> &vec,
163
				    Mesh *mesh,
164
				    string filename,
165
				    string postfix,
166 167
				    int level,
				    bool writeAsVector)
Thomas Witkowski's avatar
Thomas Witkowski committed
168
  {
169
    ElementFileWriter efw("", mesh, vec);
170 171 172 173 174 175 176 177 178 179 180 181 182
    efw.writeVtkValues(filename, postfix, level, writeAsVector);
  }


  void ElementFileWriter::writeFile(map<int, vector<double> > &vecs,
				    Mesh *mesh,
				    string filename,
				    string postfix,
				    int level,
				    bool writeAsVector)
  {
    ElementFileWriter efw("", mesh, vecs);
    efw.writeVtkValues(filename, postfix, level, writeAsVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
183 184 185
  }
  

186
  void ElementFileWriter::writeMeshDatValues(string filename, double time)
187 188
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
189

190 191 192 193 194 195 196 197 198
    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));
199

200 201
    int dim = mesh->getDim();
    double val;
202

203
    // === Write header. ===
204 205 206 207 208 209
    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";
210 211


212
    // === Write vertex coordinates (every vertex for every element). ===
213
    file << "vertex coordinates:\n";
214
    TraverseStack stack;
215

216 217 218 219
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
220

221
    while (elInfo) {
222

223 224 225
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; ++i) {     
	for (int j = 0; j < dim; ++j) {
226
	  file << elInfo->getCoord(i)[j] << " ";
227
	}
228
	file << "\n";
229 230 231 232
      }
    
      elInfo = stack.traverseNext(elInfo);
    }
233 234


235 236 237
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
238 239
    file << "\n";
    file << "element vertices:\n";
240 241
    for (int i = 0; i < numLeaves; ++i) {
      for (int j = 0; j <= dim; ++j) {
242
	file << vertCntr << " ";
243
	++vertCntr;
244
      }
245
      file << "\n";
246 247 248
    }


249
    // === Write values. ===
250

251
    // Write values header.
252 253 254 255 256 257 258 259
    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";
260 261

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

264
    file.setf(ios::scientific,ios::floatfield);
265 266 267 268 269 270 271

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

    while (elInfo) {
      // Get element value.
272
      val = (*vec)[elInfo->getElement()->getIndex()];
273 274 275
    
      // Write value for each vertex of each element.
      for (int i = 0; i <= dim; ++i) {      
276
	file << val << "\n";
277 278 279 280
      }
    
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
281 282


283
    // Write values trailor.
284 285 286
    file << "\n";
    file << "interpolation values: " << name.c_str() << "\n\n\n";
    file << "element interpolation points: " << name.c_str() << "\n";
287
  }
288

289

290 291
  void ElementFileWriter::writeVtkValues(string fname, 
					 string postfix, 
292 293
					 int level,
					 bool writeAsVector)
294 295 296
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");

297
    TEST_EXIT((vec!=NULL || vecs!=NULL) && (vec==NULL || vecs==NULL))
298
      ("Ether vec or vecs must be given, not both and not nothing!");
299

300 301 302 303 304 305 306 307 308
#if HAVE_PARALLEL_DOMAIN_AMDIS
    string filename = 
      fname  + 
      "-p" + boost::lexical_cast<std::string>(MPI::COMM_WORLD.Get_rank()) + "-" +
      postfix;
#else
    string filename = fname + postfix;
#endif

309 310 311 312 313 314 315 316 317
    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));
318 319 320

    int dim = mesh->getDim();
    int vertices = mesh->getGeo(VERTEX);
321 322
    int nVertices = 0;
    int nElements = 0;
323 324
    Flag traverseFlag = level == -1 ? Mesh::CALL_LEAF_EL : Mesh::CALL_EL_LEVEL;

325 326 327 328 329 330
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
    while (elInfo) {
      nElements++;
      elInfo = stack.traverseNext(elInfo);
    }
331

332 333 334 335 336 337 338 339 340 341 342 343 344 345
    // === Write vertex coordinates (every vertex for every element). ===
    elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);

    std::map<DegreeOfFreedom, std::pair<int, WorldVector<double> > > dof2Idx;
    std::vector<DegreeOfFreedom> idx2Dof;
    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {
	DegreeOfFreedom idx = elInfo->getElement()->getDof(i, 0);
	if (dof2Idx.find(idx) == dof2Idx.end()) {
	  dof2Idx[idx] = std::make_pair(nVertices, elInfo->getCoord(i));
	  idx2Dof.push_back(idx);
	  nVertices++;
	}
346
      }
347
      elInfo = stack.traverseNext(elInfo);
348
    }
349

350 351 352
    file << "<?xml version=\"1.0\"?>\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "  <UnstructuredGrid>\n";
353
    file << "    <Piece NumberOfPoints=\"" << nVertices
354 355 356
	 << "\" NumberOfCells=\"" <<  nElements << "\">\n";
    file << "      <Points>\n";
    file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
357

358 359 360
    for (int i = 0; i < nVertices; i++) {
      DegreeOfFreedom dof = idx2Dof[i];
      writeCoord(file, dof2Idx[dof].second);
361
    }
362

363 364 365
    file << "        </DataArray>\n";
    file << "      </Points>\n";
    file << "      <Cells>\n";
366

367
    file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
368
    for (int i = 0; i < nElements; i++) 
369 370
      file << " " << (i + 1) * vertices << "\n";    
    file << "        </DataArray>\n";
371 372


373
    file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";
374 375 376
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
377
	file << " 3\n";
378 379
	break;
      case 3: 
380
	file << " 5\n";
381 382
	break;
      case 4:
383
	file << " 10\n";
384 385 386 387 388
	break;		    
      default:
	break;
      }	   
    }
389
    file << "        </DataArray>\n";
390

391
    file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
392 393 394 395 396 397
    elInfo = stack.traverseFirst(mesh, level, traverseFlag);
    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {
	DegreeOfFreedom dof = elInfo->getElement()->getDof(i, 0);
	file << (dof2Idx[dof].first) << " ";
398
      }
399
      file << "\n";
400
      elInfo = stack.traverseNext(elInfo);
401
    }
402 403
    file << "        </DataArray>\n";
    file << "      </Cells>\n";
404

405 406 407
    int dataLength = (vecs != NULL ? (*(vecs->begin())).second.size() : 1);
    int nComponents = (!writeAsVector || (vecs == NULL && vec != NULL) ? 1 : dataLength);
    int nDataArrays = (!writeAsVector && (vec == NULL && vecs != NULL) ? dataLength : 1);
408 409 410
    file << "      <CellData>\n";
    for (int i = 0; i < nDataArrays; i++) {
      file << "        <DataArray type=\"Float32\" Name=\"value"<<i<<"\" format=\"ascii\" NumberOfComponents=\""<<(writeAsVector ? std::max(3,nComponents) : nComponents)<<"\">\n";
411

412
      file.setf(ios::scientific,ios::floatfield);
413

414 415 416 417 418 419 420
      elInfo = stack.traverseFirst(mesh, level, traverseFlag);
      int vc = 0;
      while (elInfo) {
	// Get element value.
	int idx = elInfo->getElement()->getIndex();
	
	for (int j = 0; j < nComponents; j++) {
421
	  double val = (vec != NULL ? (*vec)[idx] : (static_cast<int>((*vecs)[idx].size())==dataLength ? (*vecs)[idx][i*nComponents+j] : 0.0));
422 423 424 425 426 427

	  // Write value for each vertex of each element.
	  if (fabs(val) < 1.e-40)
	    file << "0.0";
	  else
	    file << val;
428

429 430 431
	  if (j < nComponents-1)
	    file << " ";
	}
432
	if (writeAs3dVector && writeAsVector && vecs != NULL) {
433 434 435 436 437 438 439 440 441 442 443 444 445
	  for (int j = nComponents; j < 3; j++)
	    file << " 0.0";
	}
	file << "\n";

	vc++;
	elInfo = stack.traverseNext(elInfo);
      }

      file << "        </DataArray>\n";
    }

    file << "      </CellData>\n";
446 447 448
    file << "    </Piece>\n";
    file << "  </UnstructuredGrid>\n";
    file << "</VTKFile>\n";
449 450

#if HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
451 452 453
    if (MPI::COMM_WORLD.Get_rank() == 0) {
      vector<string> componentNames;
      componentNames.push_back("elvalue");
454
      VtkWriter::detail::writeParallelFile(fname + ".pvtu", MPI::COMM_WORLD.Get_size(),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
455 456
				   fname, ".vtu", componentNames);
    }
457
#endif
458
  }
459 460
  
} } // end namespace io, AMDiS
461 462