ElementFileWriter.cc 9.22 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 17 18
#include "ElementFileWriter.h"
#include "BasisFunction.h"
#include "Parameters.h"
#include "Traverse.h"
#include "AdaptInfo.h"

19 20
namespace AMDiS {

21
  ElementFileWriter::ElementFileWriter(std::string name_, 
22
				       Mesh *mesh_,
Thomas Witkowski's avatar
Thomas Witkowski committed
23
				       std::map<int, double> &mapvec)
24 25 26 27 28 29 30 31 32 33
    : name(name_),
      amdisMeshDatExt(".elem.mesh"),
      vtkExt(".vtu"),
      writeAMDiSFormat(0),
      writeVtkFormat(0),
      appendIndex(0),
      indexLength(5),
      indexDecimals(3),
      tsModulo(1),
      timestepNumber(-1),
34
      mesh(mesh_),
Thomas Witkowski's avatar
Thomas Witkowski committed
35
      vec(mapvec)
36
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
37 38 39 40 41 42 43 44 45 46
    if (name != "") {
      GET_PARAMETER(0, name + "->output->filename", &filename);
      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);
    }
47
  }
48

49

50
  void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
Thomas Witkowski's avatar
Thomas Witkowski committed
51 52
				     int level, Flag traverseFlag, 
				     bool (*writeElem)(ElInfo*))
53 54
  {
    FUNCNAME("ElementFileWriter::writeFiles()");
55

56 57
    timestepNumber++;
    timestepNumber %= tsModulo;
58

59
    if (timestepNumber != 0 && !force) 
60
      return;
61

62
    std::string fn = filename;
63

64 65 66 67
    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");
68
    
69 70
      char formatStr[9];
      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
71

72 73
      char timeStr[20];
      sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
74

75 76
      fn += timeStr;
    }
77

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

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

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

94

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

105
  void ElementFileWriter::writeMeshDatValues(std::string filename, double time)
106 107
  {
    FUNCNAME("ElementFileWriter::writeMeshDatValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    std::ofstream fout(filename.c_str());
109

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

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

115 116 117 118 119 120 121
    // === 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";
122 123


124 125 126
    // === Write vertex coordinates (every vertex for every element). ===
    fout << "vertex coordinates:\n";
    TraverseStack stack;
127

128 129 130 131
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
132

133
    while (elInfo) {
134

135 136 137 138 139 140 141 142 143 144
      // 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);
    }
145 146


147 148 149 150 151 152 153 154 155
    // === 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;
156 157 158 159 160
      }
      fout << "\n";
    }


161
    // === Write values. ===
162

163
    // Write values header.
164
    fout << "\n";
165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
    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
193 194


195 196 197 198
    // Write values trailor.
    fout << "\n";
    fout << "interpolation values: " << name.c_str() << "\n\n\n";
    fout << "element interpolation points: " << name.c_str() << "\n";
199 200


201 202
    fout.close();
  }
203

204

205
  void ElementFileWriter::writeVtkValues(std::string filename, int level)
206 207
  {
    FUNCNAME("ElementFileWriter::writeVtkValues()");
Thomas Witkowski's avatar
Thomas Witkowski committed
208
    std::ofstream fout(filename.c_str());
209

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

    int dim = mesh->getDim();
213
    int dimOfWorld = Global::getGeo(WORLD);
214 215
    int vertices = mesh->getGeo(VERTEX);
    int nElements = mesh->getNumberOfLeaves();
216 217 218 219 220 221 222 223 224 225 226 227
    Flag traverseFlag = level == -1 ? Mesh::CALL_LEAF_EL : Mesh::CALL_EL_LEVEL;

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

      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, level, Mesh::CALL_EL_LEVEL);      
      while (elInfo) {
	nElements++;
	elInfo = stack.traverseNext(elInfo);
      }
    }
228

229 230 231
    fout << "<?xml version=\"1.0\"?>" << std::endl;
    fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
    fout << "  <UnstructuredGrid>" << std::endl;
232
    fout << "    <Piece NumberOfPoints=\"" << (dim + 1) * nElements
233 234 235
	 << "\" NumberOfCells=\"" <<  nElements << "\">" << std::endl;
    fout << "      <Points>" << std::endl;
    fout << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
236 237 238 239


    // === Write vertex coordinates (every vertex for every element). ===
    TraverseStack stack;
240
    ElInfo *elInfo = 
241
      stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
242 243 244 245

    while (elInfo) {
      // Write coordinates of all element vertices.
      for (int i = 0; i <= dim; i++) {     
246
	for (int j = 0; j < dimOfWorld; j++) 
247
	  fout << elInfo->getCoord(i)[j] << " ";
248 249
	
	for (int j = dimOfWorld; j < 3; j++) 
250
	  fout << "0.0";
251
	
252 253 254 255 256
	fout << "\n";
      }
   
      elInfo = stack.traverseNext(elInfo);
    }
257

258 259 260
    fout << "        </DataArray>" << std::endl;
    fout << "      </Points>" << std::endl;
    fout << "      <Cells>" << std::endl;
261

262
    fout << "        <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
263 264
    for (int i = 0; i < nElements; i++) 
      fout << " " << (i + 1) * vertices << std::endl;    
265
    fout << "        </DataArray>" << std::endl;
266 267


268
    fout << "        <DataArray type=\"UInt8\" Name=\"types\">" << std::endl;
269 270 271
    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
272
	fout << " 3" << std::endl;
273 274
	break;
      case 3: 
275
	fout << " 5" << std::endl;
276 277
	break;
      case 4:
278
	fout << " 10" << std::endl;
279 280 281 282 283
	break;		    
      default:
	break;
      }	   
    }
284
    fout << "        </DataArray>" << std::endl;
285

286
    fout << "        <DataArray type=\"Int32\" Name=\"connectivity\">" << std::endl;
287 288 289 290 291 292
    int vertCntr = 0;
    for (int i = 0; i < nElements; ++i) {
      for (int j = 0; j <= dim; ++j) {
	fout << vertCntr << " ";
	++vertCntr;
      }
293
      fout << std::endl;
294
    }
295
    fout << "        </DataArray>" << std::endl;
296
    
297 298 299
    fout << "      </Cells>" << std::endl;
    fout << "      <PointData>" << std::endl;
    fout << "        <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">" << std::endl;
300 301 302

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

303
    elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
304 305 306
    int vc = 0;
    while (elInfo) {
      // Get element value.
307
      double val = vec[elInfo->getElement()->getIndex()];
308
   
309
      // Write value for each vertex of each element.
310
      for (int i = 0; i <= dim; i++) 
311
	fout << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
312
      
313 314 315
      vc++;
      elInfo = stack.traverseNext(elInfo);
    } 
316

317
    fout << "        </DataArray>" << std::endl;
318
    
319

320 321 322 323
    fout << "      </PointData>" << std::endl;
    fout << "    </Piece>" << std::endl;
    fout << "  </UnstructuredGrid>" << std::endl;
    fout << "</VTKFile>" << std::endl;
324

325 326
    fout.close();
  }
327
}