ElementFileWriter.cc 9.98 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13 14 15 16
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/convenience.hpp>
17
#include <boost/lexical_cast.hpp>
18

19
#include "io/VtkWriter.h"
20 21
#include "ElementFileWriter.h"
#include "BasisFunction.h"
22
#include "Initfile.h"
23 24 25
#include "Traverse.h"
#include "AdaptInfo.h"

26 27
namespace AMDiS {

28 29 30
  using namespace std;

  ElementFileWriter::ElementFileWriter(string name_, 
31
				       Mesh *mesh_,
32
				       map<int, double> &mapvec)
33 34 35 36 37 38 39 40 41 42
    : name(name_),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
      writeAMDiSFormat(0),
      writeVtkFormat(0),
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
43
      mesh(mesh_),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      vec(mapvec)
45
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
46
    if (name != "") {
47 48 49 50 51 52 53 54
      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->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
55
    }
56
  }
57

58

59
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
60 61
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
62 63
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
64

65 66
    timestepNumber++;
    timestepNumber %= tsModulo;
67

68
    if (timestepNumber != 0 && !force) 
69
      return;
70

71
    string fn = filename;
72

73 74 75 76
    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");
77
    
78 79
      char formatStr[9];
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
80

81 82
      char timeStr[20];
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
83

84 85
      fn += timeStr;
    }
86

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

95 96
    if (writeVtkFormat) {
      TEST_EXIT(mesh)("no mesh\n");
97

98
      writeVtkValues(fn, vtkExt);
99 100
      MSG("VTK file written to %s\n", (fn + vtkExt).c_str());		     
    }
101 102
  }

103

104
  void ElementFileWriter::writeFile(map<int, double> &vec,
105
				    Mesh *mesh,
106
				    string filename,
107
				    string postfix,
108
				    int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
109
  {
110
    ElementFileWriter efw("", mesh, vec);
111
    efw.writeVtkValues(filename, postfix, level);
Thomas Witkowski's avatar
Thomas Witkowski committed
112 113 114
  }
  

115
  void ElementFileWriter::writeMeshDatValues(string filename, double time)
116 117
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
118

119 120 121 122 123 124 125 126 127
    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));
128

129 130
    int dim = mesh->getDim();
    double val;
131

132
    // === Write header. ===
133 134 135 136 137 138
    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";
139 140


141
    // === Write vertex coordinates (every vertex for every element). ===
142
    file << "vertex coordinates:\n";
143
    TraverseStack stack;
144

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

150
    while (elInfo) {
151

152 153 154
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; ++i) {     
	for (int j = 0; j < dim; ++j) {
155
	  file << elInfo->getCoord(i)[j] << " ";
156
	}
157
	file << "\n";
158 159 160 161
      }
    
      elInfo = stack.traverseNext(elInfo);
    }
162 163


164 165 166
    // === Write elements. ===
    int numLeaves = mesh->getNumberOfLeaves();
    int vertCntr = 0;
167 168
    file << "\n";
    file << "element vertices:\n";
169 170
    for (int i = 0; i < numLeaves; ++i) {
      for (int j = 0; j <= dim; ++j) {
171
	file << vertCntr << " ";
172
	++vertCntr;
173
      }
174
      file << "\n";
175 176 177
    }


178
    // === Write values. ===
179

180
    // Write values header.
181 182 183 184 185 186 187 188
    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";
189 190

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

193
    file.setf(ios::scientific,ios::floatfield);
194 195 196 197 198 199 200 201 202 203 204

    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) {      
205
	file << val << "\n";
206 207 208 209
      }
    
      elInfo = stack.traverseNext(elInfo);
    }  // end of: mesh traverse
210 211


212
    // Write values trailor.
213 214 215
    file << "\n";
    file << "interpolation values: " << name.c_str() << "\n\n\n";
    file << "element interpolation points: " << name.c_str() << "\n";
216
  }
217

218

219 220 221
  void ElementFileWriter::writeVtkValues(string fname, 
					 string postfix, 
					 int level)
222 223 224
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");

225 226 227 228 229 230 231 232 233
#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

234 235 236 237 238 239 240 241 242
    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));
243 244

    int dim = mesh->getDim();
245
    int dimOfWorld = Global::getGeo(WORLD);
246 247
    int vertices = mesh->getGeo(VERTEX);
    int nElements = mesh->getNumberOfLeaves();
248 249 250 251 252 253
    Flag traverseFlag = level == -1 ? Mesh::CALL_LEAF_EL : Mesh::CALL_EL_LEVEL;

    if (level != -1) {
      nElements = 0;

      TraverseStack stack;
254
      ElInfo *elInfo = stack.traverseFirst(mesh, level, Mesh::CALL_EL_LEVEL);
255 256 257 258 259
      while (elInfo) {
	nElements++;
	elInfo = stack.traverseNext(elInfo);
      }
    }
260

261 262 263 264 265 266 267
    file << "<?xml version=\"1.0\"?>\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "  <UnstructuredGrid>\n";
    file << "    <Piece NumberOfPoints=\"" << (dim + 1) * nElements
	 << "\" NumberOfCells=\"" <<  nElements << "\">\n";
    file << "      <Points>\n";
    file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
268 269 270 271


    // === Write vertex coordinates (every vertex for every element). ===
    TraverseStack stack;
272
    ElInfo *elInfo = 
273
      stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
274 275 276 277

    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {     
278
	for (int j = 0; j < dimOfWorld; j++) 
279
	  file << elInfo->getCoord(i)[j] << " ";
280 281
	
	for (int j = dimOfWorld; j < 3; j++) 
282
	  file << "0.0";
283
	
284
	file << "\n";
285 286 287 288
      }
   
      elInfo = stack.traverseNext(elInfo);
    }
289

290 291 292
    file << "        </DataArray>\n";
    file << "      </Points>\n";
    file << "      <Cells>\n";
293

294
    file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
295
    for (int i = 0; i < nElements; i++) 
296 297
      file << " " << (i + 1) * vertices << "\n";    
    file << "        </DataArray>\n";
298 299


300
    file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";
301 302 303
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
304
	file << " 3\n";
305 306
	break;
      case 3: 
307
	file << " 5\n";
308 309
	break;
      case 4:
310
	file << " 10\n";
311 312 313 314 315
	break;		    
      default:
	break;
      }	   
    }
316
    file << "        </DataArray>\n";
317

318
    file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
319 320 321
    int vertCntr = 0;
    for (int i = 0; i < nElements; ++i) {
      for (int j = 0; j <= dim; ++j) {
322
	file << vertCntr << " ";
323 324
	++vertCntr;
      }
325
      file << "\n";
326
    }
327
    file << "        </DataArray>\n";
328
    
329 330 331
    file << "      </Cells>\n";
    file << "      <PointData>\n";
    file << "        <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">\n";
332

333
    file.setf(ios::scientific,ios::floatfield);
334

335
    elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
336 337 338
    int vc = 0;
    while (elInfo) {
      // Get element value.
339
      double val = vec[elInfo->getElement()->getIndex()];
340
   
341
      // Write value for each vertex of each element.
342
      for (int i = 0; i <= dim; i++) 
343
	file << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
344
      
345 346 347
      vc++;
      elInfo = stack.traverseNext(elInfo);
    } 
348

349
    file << "        </DataArray>\n";
350
    
351

352 353 354 355
    file << "      </PointData>\n";
    file << "    </Piece>\n";
    file << "  </UnstructuredGrid>\n";
    file << "</VTKFile>\n";
356 357 358 359 360 361

#if HAVE_PARALLEL_DOMAIN_AMDIS
    if (MPI::COMM_WORLD.Get_rank() == 0)
      VtkWriter::writeParallelFile(fname + ".pvtu", MPI::COMM_WORLD.Get_size(),
				   fname, ".vtu", 1);
#endif
362
  }
363
}