Arh3Writer.cc 26.2 KB
Newer Older
Siqi Ling's avatar
Siqi Ling committed
1
2
3
#include <fstream>
#include <stdint.h>
#include <iostream>
4
#include <streambuf>
Siqi Ling's avatar
Siqi Ling committed
5
6
7
8
9
10
11
12

#include "Arh3Writer.h"
#include "Mesh.h"
#include "MeshStructure.h"
#include "Traverse.h"
#include "DOFVector.h"
#include "../Arh3Reader.h"

13
#include <boost/filesystem.hpp>
Siqi Ling's avatar
Siqi Ling committed
14
#include <boost/iostreams/filtering_stream.hpp>
Siqi Ling's avatar
Siqi Ling committed
15
#include <boost/iostreams/filtering_streambuf.hpp>
Siqi Ling's avatar
Siqi Ling committed
16
#include <boost/iostreams/device/file_descriptor.hpp>
Siqi Ling's avatar
Siqi Ling committed
17
#include <boost/iostreams/copy.hpp>
18
#ifdef AMDIS_HAS_COMPRESSION
Siqi Ling's avatar
Siqi Ling committed
19
20
21
22
23
#include <boost/iostreams/filter/zlib.hpp>
#include <boost/iostreams/filter/bzip2.hpp>
#endif

namespace AMDiS { namespace io {
24

Siqi Ling's avatar
Siqi Ling committed
25
26
27
28
29
30
  using namespace std;

  namespace Arh3Writer
  {
    namespace detail
    {
31
32
      void write(string filename,
  			DOFVector<double>* vec0,
Siqi Ling's avatar
Siqi Ling committed
33
34
35
  			DOFVector<double>* vec1,
  			DOFVector<double>* vec2,
			bool writeParallel,
36
			Cpsformat::Value cps,
Siqi Ling's avatar
Siqi Ling committed
37
			string dataformat,
38
			Macroformat::Value writeMacro)
Siqi Ling's avatar
Siqi Ling committed
39
40
41
42
43
44
45
46
47
      {
        vector<DOFVector<double>*> vecs(0);
        if (vec0 != NULL)
          vecs.push_back(vec0);
        if (vec1 != NULL)
          vecs.push_back(vec1);
        if (vec2 != NULL)
          vecs.push_back(vec2);

48
        write(filename, NULL, vecs, writeParallel, cps, dataformat, writeMacro);
Siqi Ling's avatar
Siqi Ling committed
49
      }
50
51

      void updateAnimationFile(AdaptInfo *adaptInfo,
Siqi Ling's avatar
Siqi Ling committed
52
53
54
55
56
57
58
59
60
61
62
63
			       std::vector<double>* arhAnimationFrames,
			       bool createSubDir,
			       int indexLength,
			       int indexDecimals,
			       std::string animationFilename)
      {
	FUNCNAME("updateAnimationFile()");

	arhAnimationFrames->push_back(adaptInfo->getTime());

	boost::iostreams::filtering_ostream file;
	{
64
	  ofstream swapfile(animationFilename.c_str(),
Siqi Ling's avatar
Siqi Ling committed
65
66
67
68
69
			    ios::out | ios::trunc | ios::binary);
	  TEST_EXIT(swapfile.is_open())
	    ("Cannot open file %s for writing!\n", animationFilename.c_str());
	  swapfile.close();
	}
70
	file.push(boost::iostreams::file_descriptor_sink(animationFilename,
Siqi Ling's avatar
Siqi Ling committed
71
							  ios::trunc | ios::binary));
72

Siqi Ling's avatar
Siqi Ling committed
73
74
75
76
77
78
79
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
	string typeId = "tparh";
#else
	string typeId = "tarh";
#endif
	string baseDir = createSubDir ?  "./data/" : "./";
	uint32_t baseDirLen = baseDir.length();
80
	int major = AMDiS::io::Arh3Reader::TARH_MAJOR;
Siqi Ling's avatar
Siqi Ling committed
81
82
83
84
85
86
87
88
        int minor = AMDiS::io::Arh3Reader::TARH_MINOR;
	typeId += '_' + boost::lexical_cast<string>(major) + '.' + boost::lexical_cast<string>(minor);
	if (typeId.length() <= 16) {	// 16 is the Id size
	  string rest(16 - typeId.length(), ' ');
	  typeId.append(rest);
	}
	uint32_t indexLength_ = indexLength;
	uint32_t indexDecimals_ = indexDecimals;
89

Siqi Ling's avatar
Siqi Ling committed
90
91
92
93
94
	file.write(typeId.c_str(), 16);
	file.write(reinterpret_cast<char*>(&baseDirLen), 4);
	file.write(baseDir.c_str(), baseDirLen);
	file.write(reinterpret_cast<char*>(&indexLength_), 4);
	file.write(reinterpret_cast<char*>(&indexDecimals_), 4);
95

Siqi Ling's avatar
Siqi Ling committed
96
97
98
	for (size_t i = 0; i < arhAnimationFrames->size(); i++)
	  file.write(reinterpret_cast<char*>(&arhAnimationFrames[i]), 8);
      }
99
100

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
101
      void writeParallelFile(string filename, Mesh* mesh, bool createSubDir, Macroformat::Value writeMacro)
Siqi Ling's avatar
Siqi Ling committed
102
103
104
      {
	ofstream file;
        file.open(filename.c_str(), ios::out | ios::binary | ios::trunc);
105

106
	string typeId = "parh", macroFilename = "", perFilename = "";
107
	string baseDir = createSubDir ?  "./data" : ".";
108
109
	string macroFile = "";
	uint32_t macroFile_nl = 0;
110
111
	string macroData = "", periodicData = "";
	uint32_t baseDirLen = baseDir.length();
Siqi Ling's avatar
Siqi Ling committed
112
	Parameters::get(mesh->getName() + "->macro file name", macroFilename);
113
	Parameters::get(mesh->getName() + "->periodic file", perFilename);
114

115
	bool hasPeriodic = perFilename.length();
116
	int major = AMDiS::io::Arh3Reader::PARH_MAJOR;
Siqi Ling's avatar
Siqi Ling committed
117
        int minor = AMDiS::io::Arh3Reader::PARH_MINOR;
Siqi Ling's avatar
Siqi Ling committed
118
	uint32_t nFiles = MPI::COMM_WORLD.Get_size();
119
	map<int, int> partitionMap =
Siqi Ling's avatar
Siqi Ling committed
120
	  Parallel::MeshDistributor::globalMeshDistributor->getPartitionMap();
121
	uint32_t nMacros = partitionMap.size();
122

123
124
125
126
127
128
129
	typeId += '_' + boost::lexical_cast<string>(major) + '.' + boost::lexical_cast<string>(minor);
	if (typeId.length() <= 16) {	// 16 is the Id size
	  string rest(16 - typeId.length(), ' ');
	  typeId.append(rest);
	}
	else
	  ERROR_EXIT("Should not happen.\n");
130

131
132
133
	file.write(typeId.c_str(), 16);
	file.write(reinterpret_cast<char*>(&baseDirLen), 4);
	file.write(baseDir.c_str(), baseDirLen);
Siqi Ling's avatar
Siqi Ling committed
134
	file.write(reinterpret_cast<char*>(&nFiles), 4);
135

136
	if (writeMacro != Macroformat::NONE && macroFilename.length()) {
137

138
	  if (writeMacro == Macroformat::PT_MACROFILE) {
139
	    macroFile = !hasPeriodic ? macroFilename :
140
	      macroFilename + ';' + perFilename;
141

142
	    macroFile_nl = macroFile.length();
143
	  } else if (writeMacro == Macroformat::SELF) {
144
145
	    macroFile_nl = 13;
	    readFileToString(macroFilename, macroData);
146

147
148
149
150
151
	    if (hasPeriodic) {
	      macroFile_nl = 27;
	      readFileToString(perFilename, periodicData);
	    }
	    macroFile.resize(macroFile_nl, ' ');
152
153
	  }
	}
Siqi Ling's avatar
Siqi Ling committed
154
	file.write(reinterpret_cast<char*>(&macroFile_nl), 4);
155
	file.write(macroFile.c_str(), macroFile_nl);
156
	file.write(reinterpret_cast<char*>(&nMacros), 4);
157

Siqi Ling's avatar
Siqi Ling committed
158
159
160
161
162
163
	map<int, int>::const_iterator it = partitionMap.begin();
	uint32_t rank = 0;
	for (;it != partitionMap.end(); it++) {
	  rank = it->second;
	  file.write(reinterpret_cast<char*>(&rank), 4);
	}
164

165
	// write macro and periodic file
166
	if (writeMacro == Macroformat::SELF && macroFilename.length()) {
167
168
169
	  file.seekp(0, ios_base::end);
	  long macroPos = file.tellp(), perPos = 0;
	  file.write(macroData.c_str(), macroData.length());
170

171
	  if (hasPeriodic) {
172
173
174
	    perPos = file.tellp();
	    file.write(periodicData.c_str(), periodicData.length());
	  }
175

176
177
178
179
180
181
182
	  // update macroFile_nl
	  int offset = 16 +		//typeId
		       4 +		//baseDirLen
		       baseDirLen +	//baseDir
		       4 +		//nFiles
		       4;		//macroFile_nl
	  file.seekp(offset);
183
184
185
186
187
188
	  file.write("this:", 5);
	  file.write(reinterpret_cast<char*>(&macroPos), 8);
	  if (hasPeriodic) {
	    file.write(";this:", 6);
	    file.write(reinterpret_cast<char*>(&perPos), 8);
	  }
189
	}
Siqi Ling's avatar
Siqi Ling committed
190
191
192
      }
#endif

193
194
      void write(string filename,
		 Mesh* mesh,
Siqi Ling's avatar
Siqi Ling committed
195
196
		 vector<DOFVector<double>*> vecs,
		 bool writeParallel,
197
		 Cpsformat::Value cps,
Siqi Ling's avatar
Siqi Ling committed
198
		 string dataformat,
199
		 Macroformat::Value writeMacro)
Siqi Ling's avatar
Siqi Ling committed
200
201
      {
	FUNCNAME("Arh3Writer::detail::write()");
202

203
204
	int sPos = filename.find(".arh");
	TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n");
205

Siqi Ling's avatar
Siqi Ling committed
206
207
	map<string,Valformat>::const_iterator it = dataformatMap.find(dataformat);
	TEST_EXIT(it != dataformatMap.end())("Wrong data format.\n");
208

209
210
211
212
213
	std::set<Mesh*> meshset;
	std::set<string> nameset;
	if (mesh)
	  meshset.insert(mesh);
	for (size_t i = 0; i < vecs.size(); i++) {
Siqi Ling's avatar
Siqi Ling committed
214
	  TEST_EXIT(vecs[i] != NULL)("Vecs[%i] is NULL. Please check.\n", i);
215
216
217
	  meshset.insert(vecs[i]->getFeSpace()->getMesh());
	  nameset.insert(vecs[i]->getName());
	}
218

219
220
221
	if (meshset.size() == 0) {
	  WARNING("There is nothing to be writen.\n");
	  return;
Siqi Ling's avatar
Siqi Ling committed
222
	}
223

224
225
	TEST_EXIT(nameset.size() == vecs.size())
	  ("DOFVectors in vecs cannot have idential name. Please check.\n");
226
227
228
229

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
	if (writeParallel) {
	  using boost::lexical_cast;
230
	  using boost::filesystem::path;
231

232
	  Mesh* mesh_ = mesh ? mesh : vecs[0]->getFeSpace()->getMesh();
233

234
235
	  bool arh3CreateSubDur = false;
	  Parameters::get("arh3->create subdirectory", arh3CreateSubDur);
236

237
238
	  string parhname = filename.substr(0, sPos);
	  string sarhname = parhname;
239

240
241
242
243
244
245
	  if (arh3CreateSubDur) {
	    path vtu_path = parhname;
	    path data_basedir("data");
	    path vtu_filename = vtu_path.filename();
	    vtu_path.remove_filename() /= data_basedir;
	    try {
246
247
	      create_directory(vtu_path);
	      vtu_path /= vtu_filename;
248
249
250
	      sarhname = vtu_path.string();
	    } catch (...) {}
	  }
251

252
	  if (MPI::COMM_WORLD.Get_rank() == 0)
253
	    writeParallelFile(parhname + ".parh", mesh_, arh3CreateSubDur);
254

255
	  filename = sarhname + "_p" + lexical_cast<string>(MPI::COMM_WORLD.Get_rank()) + ".arh";
256
	}
257
258
#endif

259
	bool multiMesh = meshset.size() > 1;
260

Siqi Ling's avatar
Siqi Ling committed
261
	//if mesh exists, the meshes in vecs should be the same.
262
263
	if (mesh) {
	  for (size_t i = 0; i < vecs.size(); i++)
Siqi Ling's avatar
Siqi Ling committed
264
	    TEST_EXIT(mesh == vecs[i]->getFeSpace()->getMesh())
265
	      ("The mesh of DOFVector %i in vecs is not equal to the second parameter.\n", i);
266
	  writeAux(filename, mesh, vecs, writeParallel, cps, dataformat, writeMacro);
267
268
	} else {
	  if (!multiMesh)
269
	    writeAux(filename, vecs[0]->getFeSpace()->getMesh(), vecs, writeParallel, cps, dataformat, writeMacro);
270
271
272
273
	  else {
	    vector<bool> visited(vecs.size(), false);
	    vector<DOFVector<double>*> splitedVecs(0);
	    Mesh* tmpMesh = NULL;
274

275
276
277
278
279
280
281
282
283
284
285
286
	    for(size_t i = 0; i < vecs.size(); i++) {
	      if(!visited[i]) {
		splitedVecs.clear();
		splitedVecs.push_back(vecs[i]);
		visited[i] = true;
		tmpMesh = vecs[i]->getFeSpace()->getMesh();
		for(size_t j = i + 1; j < vecs.size(); j++)
		  if(vecs[j]->getFeSpace()->getMesh() == tmpMesh)
		  {
		    splitedVecs.push_back(vecs[j]);
		    visited[j] = true;
		  }
Siqi Ling's avatar
Siqi Ling committed
287
288
	      }
	      string newfilename = filename;
289
290
291
292
293
	      if(filename.length() > 4 && filename.substr(filename.length()- 4, filename.length()) == ".arh")
		newfilename = filename.substr(0, filename.length() - 4) +
		  "." + tmpMesh->getName() + filename.substr(filename.length()-4 , filename.length());
	      else
		newfilename = filename + "." + tmpMesh->getName() + ".arh";
294
295
296
297
298
299

	      writeAux(newfilename,
		       splitedVecs[0]->getFeSpace()->getMesh(),
		       splitedVecs,
		       writeParallel,
		       cps,
300
		       dataformat,
301
		       writeMacro);
Siqi Ling's avatar
Siqi Ling committed
302
303
	    }
	  }
304
        }
Siqi Ling's avatar
Siqi Ling committed
305
      }
306

Siqi Ling's avatar
Siqi Ling committed
307
308
309
310
      int writeHeader(ofstream& file,
                               Mesh *mesh,
                               vector<DOFVector<double>*> vecs,
                               map<const FiniteElemSpace*, vector<int> >& feSpaces,
311
312
313
			       Cpsformat::Value cps,
			       string dataformat,
			       Macroformat::Value writeMacro)
Siqi Ling's avatar
Siqi Ling committed
314
315
      {
        FUNCNAME("Arh3Writer::detail::writeHeader()");
316
// 	int nbits = boost::lexical_cast<int>(dataformat.substr(2, 2));
Siqi Ling's avatar
Siqi Ling committed
317
        TEST_EXIT(file.is_open())("the file is not open. should not happen.\n");
318

Siqi Ling's avatar
Siqi Ling committed
319
	map<const FiniteElemSpace*, string> AFEDfileName;
320

Siqi Ling's avatar
Siqi Ling committed
321
        uint32_t valueNamesLen = 0, feSpacesNamesLen = 0;
Siqi Ling's avatar
Siqi Ling committed
322
323
        for (size_t i = 0; i < vecs.size(); i++)
          valueNamesLen += vecs[i]->getName().length();
324

Siqi Ling's avatar
Siqi Ling committed
325
	// TODO AFEDfileName
Siqi Ling's avatar
Siqi Ling committed
326
	map<const FiniteElemSpace*, vector<int> >::iterator feSpaceIt;
Siqi Ling's avatar
Siqi Ling committed
327
	for (feSpaceIt = feSpaces.begin(); feSpaceIt != feSpaces.end(); feSpaceIt++) {
Siqi Ling's avatar
Siqi Ling committed
328
	  AFEDfileName.insert(make_pair(feSpaceIt->first, string()));
Siqi Ling's avatar
Siqi Ling committed
329
330
	  feSpacesNamesLen += 0;
	}
331

Siqi Ling's avatar
Siqi Ling committed
332
333
334
        uint32_t nValueVectors = vecs.size();
        uint32_t nFeSpaces = feSpaces.size();
	uint32_t nMacroElements = mesh->getNumberOfMacros();
335

Siqi Ling's avatar
Siqi Ling committed
336
337
338
339
340
341
342
	uint32_t nMacro = 0;
	TraverseStack st;
	ElInfo *el = st.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
	while (el) {
	  nMacro++;
	  el = st.traverseNext(el);
	}
343

344
345
346
347
348
	string macroFilename = "", perFilename = "";
	string macroFile = "";
	uint32_t macroFile_nl = 0;
	Parameters::get(mesh->getName() + "->macro file name", macroFilename);
	Parameters::get(mesh->getName() + "->periodic file", perFilename);
349

350
	if (writeMacro != Macroformat::NONE && macroFilename.length()) {
351

352
	  if (writeMacro == Macroformat::PT_MACROFILE) {
353
	    macroFile = !perFilename.length() ? macroFilename :
354
	      macroFilename + ';' + perFilename;
355

356
	    macroFile_nl = macroFile.length();
357
	  } else if (writeMacro == Macroformat::SELF) {
358
359
360
361
362
363
	    macroFile_nl = perFilename.length() ? 27 : 13;
	    macroFile.resize(macroFile_nl, ' ');
	  }
	}

	uint32_t dow = mesh->getGeo(WORLD);
Siqi Ling's avatar
Siqi Ling committed
364
        uint32_t dim = mesh->getDim();
365
        uint32_t headerLen = 34 +                         //fixed part of header
Siqi Ling's avatar
Siqi Ling committed
366
			     nMacroElements * 12 + 12 +   //macroElemnts table
Siqi Ling's avatar
Siqi Ling committed
367
			     feSpacesNamesLen +		  //feSpaces table
Siqi Ling's avatar
Siqi Ling committed
368
369
			     nFeSpaces * 20 + 	          //feSpaces table
			     valueNamesLen + 	          //value vector table
370
			     nValueVectors * 12 +         //also value vector table
371
372
			     4 +			  //macroFile_nl
			     macroFile_nl;		  //macroFile
Siqi Ling's avatar
Siqi Ling committed
373
        string typeId = "sarh";
374
#if !AMDIS_HAS_COMPRESSION
375
	cps = Cpsformat::NONE;
Siqi Ling's avatar
Siqi Ling committed
376
#endif
377
        uint8_t *major = const_cast<uint8_t*>(&(AMDiS::io::Arh3Reader::MAJOR));
Siqi Ling's avatar
Siqi Ling committed
378
379
380
        uint8_t *minor = const_cast<uint8_t*>(&(AMDiS::io::Arh3Reader::MINOR));
	int cpsflag = static_cast<int>(cps);
	uint32_t minus1 = -1;
381

Siqi Ling's avatar
Siqi Ling committed
382
383
384
385
386
387
388
389
390
391
392
393
394
	//fixed header
        file.write(typeId.c_str(), 4);
        file.write(reinterpret_cast<char*>(major), 1);
        file.write(reinterpret_cast<char*>(minor), 1);
        file.write(reinterpret_cast<char*>(&headerLen), 4);
        file.write(reinterpret_cast<char*>(&dow), 4);
        file.write(reinterpret_cast<char*>(&dim), 4);
        file.write(reinterpret_cast<char*>(&nFeSpaces), 4);
        file.write(reinterpret_cast<char*>(&nValueVectors), 4);
        file.write(reinterpret_cast<char*>(&nMacroElements), 4);
	file.write(reinterpret_cast<char*>(&cpsflag), 4);
	//macro table
	deque<MacroElement*>::const_iterator macroIter = mesh->firstMacroElement();
395

Siqi Ling's avatar
Siqi Ling committed
396
397
398
399
400
401
402
403
404
405
406
	while(macroIter != mesh->endOfMacroElements())
	{
	  uint32_t macroIndex = (*macroIter)->getIndex();
	  file.write(reinterpret_cast<char*>(&macroIndex), 4);
	  file.write(reinterpret_cast<char*>(&minus1), 4);
	  file.write(reinterpret_cast<char*>(&minus1), 4);
	  macroIter++;
	}
	file.write(reinterpret_cast<char*>(&minus1), 4);
	file.write(reinterpret_cast<char*>(&minus1), 4);
	file.write(reinterpret_cast<char*>(&minus1), 4);
407

Siqi Ling's avatar
Siqi Ling committed
408
409
410
411
	vector<int> feSpaceNumOfVecs(vecs.size());
	uint32_t posDOFs = 0, nameStrLen = 0;
	string nameStr("");
	size_t i = 0;
412

Siqi Ling's avatar
Siqi Ling committed
413
414
415
416
417
418
419
420
421
	//feSpace table
	for(feSpaceIt = feSpaces.begin(); feSpaceIt != feSpaces.end(); feSpaceIt++, i++)
	{
	  nameStr = AFEDfileName[feSpaceIt->first];
	  nameStrLen = nameStr.length();
	  file.write(reinterpret_cast<char*>(&nameStrLen), 4);
	  file.write(nameStr.c_str(), nameStrLen);
	  DimVec<int>* nDOF = feSpaceIt->first->getBasisFcts()->getNumberOfDofs();
	  //
422
	  for(int j = 1; j < nDOF->getSize(); j++) {
Siqi Ling's avatar
Siqi Ling committed
423
424
425
            posDOFs = (*nDOF)[j];
            file.write(reinterpret_cast<char*>(&posDOFs), 4);
          }
426
          for(int j = nDOF->getSize(); j < 4 ; j++) {
Siqi Ling's avatar
Siqi Ling committed
427
428
429
430
431
432
433
434
435
            posDOFs = 0;
            file.write(reinterpret_cast<char*>(&posDOFs), 4);
          }
          posDOFs = (*nDOF)[0];
	  file.write(reinterpret_cast<char*>(&posDOFs), 4);
	  //
          for(size_t j = 0; j < feSpaceIt->second.size(); j++)
	    feSpaceNumOfVecs[feSpaceIt->second[j]] = i;
	}
436

Siqi Ling's avatar
Siqi Ling committed
437
438
439
440
441
442
443
444
445
446
	//vector table
        for(i = 0; i < vecs.size(); i++)
        {
          nameStr = vecs[i]->getName();
          nameStrLen = nameStr.length();
          file.write(reinterpret_cast<char*>(&nameStrLen), 4);
          file.write(nameStr.c_str(), nameStrLen);
	  file.write(reinterpret_cast<char*>(&feSpaceNumOfVecs[i]), 4);
	  file.write(dataformat.c_str(), 4);
        }
447

448
        file.write(reinterpret_cast<char*>(&macroFile_nl), 4);
449
	file.write(macroFile.c_str(), macroFile_nl);
Siqi Ling's avatar
Siqi Ling committed
450
451
        return headerLen;
      }
452

Siqi Ling's avatar
Siqi Ling committed
453
454
455
      void writeAux(string filename, Mesh *mesh,
                        vector<DOFVector<double>*> vecs,
                        bool writeParallel,
456
			Cpsformat::Value cps,
457
			string dataformat,
458
			Macroformat::Value writeMacro)
Siqi Ling's avatar
Siqi Ling committed
459
460
      {
        FUNCNAME("Arh3Writer::detail::writeAux()");
461

462
	TEST_EXIT(mesh)("empty mesh.\n");
463

Siqi Ling's avatar
Siqi Ling committed
464
465
466
467
468
469
470
        //initialization
        ofstream file;
        file.open(filename.c_str(), ios::out | ios::binary | ios::trunc);

        map<const FiniteElemSpace*, vector<int> > sortedFeSpaces;
	map<const FiniteElemSpace*, vector<int> >::iterator feSpaceIt;
        vector<pair<int, int> > macroSize; // (uncompressed size, compressed size)
471
472

        DegreeOfFreedom globalDof;
Siqi Ling's avatar
Siqi Ling committed
473
        size_t i = 0, j = 0;
474

Siqi Ling's avatar
Siqi Ling committed
475
476
477
478
479
480
        for(i = 0; i < vecs.size(); i++)
	{
          sortedFeSpaces[vecs[i]->getFeSpace()].push_back(i);
	}
	vector<std::set<DegreeOfFreedom> > visited(sortedFeSpaces.size());
	pair<std::set<DegreeOfFreedom>::iterator,bool> ret;
481
        //file header information
482
        int headerLen = writeHeader(file, mesh, vecs, sortedFeSpaces, cps, dataformat, writeMacro);
483

Siqi Ling's avatar
Siqi Ling committed
484
485
486
487
        //macro elements information
        MeshStructure elementStructure;
        vector<vector<double> > values(vecs.size());
        int32_t macroElIndex = -1;
488

Siqi Ling's avatar
Siqi Ling committed
489
490
491
492
493
494
495
496
497
        TraverseStack stack;
        ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
        while (elInfo) {
          if (elInfo->getLevel() == 0) {
            if (macroElIndex != -1) {
              elementStructure.commit();
	      macroSize.push_back(writeMacroElement(file, elementStructure, values, sortedFeSpaces, cps, dataformat));
            }
	    elementStructure.clear();
498

Siqi Ling's avatar
Siqi Ling committed
499
	    macroElIndex = elInfo->getElement()->getIndex();
500

Siqi Ling's avatar
Siqi Ling committed
501
502
	    for (j = 0; j < vecs.size(); j++)
	      values[j].clear();
503

Siqi Ling's avatar
Siqi Ling committed
504
505
506
507
            for (j = 0; j < sortedFeSpaces.size(); j++)
              visited[j].clear();
          }
          elementStructure.insertElement(elInfo->getElement()->isLeaf());
508
509

          if (elInfo->getElement()->isLeaf())
Siqi Ling's avatar
Siqi Ling committed
510
511
512
513
514
          {
            int valuePos = 0;
            for(feSpaceIt = sortedFeSpaces.begin(), i = 0; feSpaceIt != sortedFeSpaces.end(); feSpaceIt++, i++)
            {
              DOFAdmin* admin = feSpaceIt->first->getAdmin();
515
516
              TEST_EXIT(admin)("the DOFAdmin of DOFVector is null, this should not happen.\n");

Siqi Ling's avatar
Siqi Ling committed
517
518
              int n0, nd, nd0;

519
              if ((nd = admin->getNumberOfDofs(VERTEX)))  {
Siqi Ling's avatar
Siqi Ling committed
520
                int vertices = mesh->getGeo(VERTEX);
521
522
                nd0 = admin->getNumberOfPreDofs(VERTEX);
                n0 = mesh->getNode(VERTEX); //
Siqi Ling's avatar
Siqi Ling committed
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
                for (int n = 0; n < vertices; n++)
                  for(int d = 0; d < nd; d++)
                  {
                    globalDof = elInfo->getElement()->getDof(n0 + n, nd0 + d);
                    ret = visited[i].insert(globalDof);
                    if(ret.second)
                    {
                      for(j = 0 ; j < feSpaceIt->second.size(); j++)
                      {
                        values[valuePos + j].push_back((*vecs[feSpaceIt->second[j]])[globalDof]);
                      }
                    }
                  }
              }
              if (mesh->getDim() > 1) {
538
539
540
541
                if ((nd = admin->getNumberOfDofs(EDGE)))  {
                  int edges = mesh->getGeo(EDGE);
                  nd0 = admin->getNumberOfPreDofs(EDGE);
                  n0 = mesh->getNode(EDGE);
Siqi Ling's avatar
Siqi Ling committed
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
                  for (int n = 0; n < edges; n++)
                  for(int d = 0; d < nd; d++)
                  {
                    globalDof = elInfo->getElement()->getDof(n0 + n, nd0 + d);
                    ret = visited[i].insert(globalDof);
                    if(ret.second)
                    {
                      for(j = 0 ; j < feSpaceIt->second.size(); j++)
                      {
                        values[valuePos + j].push_back((*vecs[feSpaceIt->second[j]])[globalDof]);
                      }
                    }
                  }
                }
              }
              if (mesh->getDim() == 3) {
558
                if ((nd = admin->getNumberOfDofs(FACE)))  {
Siqi Ling's avatar
Siqi Ling committed
559
                  int faces = mesh->getGeo(FACE);
560
561
                  nd0 = admin->getNumberOfPreDofs(FACE);
                  n0 = mesh->getNode(FACE);
Siqi Ling's avatar
Siqi Ling committed
562
563
564
565
566
567
568
569
570
571
572
573
574
575
                  for (int n = 0; n < faces; n++)
                  for(int d = 0; d < nd; d++)
                  {
                    globalDof = elInfo->getElement()->getDof(n0 + n, nd0 + d);
                    ret = visited[i].insert(globalDof);
                    if(ret.second)
                    {
                      for(j = 0 ; j < feSpaceIt->second.size(); j++)
                      {
                        values[valuePos + j].push_back((*vecs[feSpaceIt->second[j]])[globalDof]);
                      }
                    }
                  }
                }
576
577
578
579
              }
              if ((nd = admin->getNumberOfDofs(CENTER)))  {
                nd0 = admin->getNumberOfPreDofs(CENTER);
                n0 = mesh->getNode(CENTER);
Siqi Ling's avatar
Siqi Ling committed
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
                for(int d = 0; d < nd; d++)
                {
                  globalDof = elInfo->getElement()->getDof(n0, nd0 + d);
                  ret = visited[i].insert(globalDof);
                    if(ret.second)
                    {
                      for(j = 0 ; j < feSpaceIt->second.size(); j++)
                      {
                        values[valuePos + j].push_back((*vecs[feSpaceIt->second[j]])[globalDof]);
                      }
                    }
                }
              }
              valuePos += feSpaceIt->second.size();
            }//feSpace loop
          }//isLeaf

          elInfo = stack.traverseNext(elInfo);
        }
        // And write the last macro element to file.
        TEST_EXIT_DBG(macroElIndex != -1)("Should not happen!\n");
        elementStructure.commit();
	macroSize.push_back(writeMacroElement(file, elementStructure, values, sortedFeSpaces, cps, dataformat));
	TEST_EXIT(macroSize.size() == (unsigned)mesh->getNumberOfMacros())("Should not happen.\n");
	//reset the macro positions in file
	setMacrosPos(file, headerLen, macroSize);
606

607
	if (writeMacro == Macroformat::SELF)
608
	  setMacroFile(file, headerLen, mesh, writeMacro);
609

Siqi Ling's avatar
Siqi Ling committed
610
611
612
        file.close();
        MSG("ARH file written to: %s\n", filename.c_str());
      }
613

Siqi Ling's avatar
Siqi Ling committed
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
      void setMacrosPos(ofstream& file, int headerLen,
	                    vector<pair<int, int> >& macroSize)
      {
	file.seekp(34);
	long pos = 0;
	int startPos = headerLen;
	for(size_t i = 0; i < macroSize.size(); i++)
	{
	  pos = file.tellp();
	  file.seekp(pos + 4);
	  file.write(reinterpret_cast<char*>(&startPos), 4);
	  file.write(reinterpret_cast<char*>(&macroSize[i].first), 4);
	  startPos += macroSize[i].second;
	}
	pos = file.tellp();
	file.seekp(pos + 4);
	file.write(reinterpret_cast<char*>(&startPos), 4);
      }
632

633
      void setMacroFile(std::ofstream& file, int headerLen, Mesh* mesh, Macroformat::Value writeMacro)
634
      {
635
	FUNCNAME("setMacroFile()");
636

637
638
	TEST_EXIT(writeMacro == Macroformat::SELF)
	  ("This function should only be called under Macroformat::SELF.\n");
639

640
641
	string macroFilename = "", perFilename = "";
	string macroData = "", periodicData = "";
642

643
644
	Parameters::get(mesh->getName() + "->macro file name", macroFilename);
	Parameters::get(mesh->getName() + "->periodic file", perFilename);
645
	bool hasPeriodic = perFilename.length();
646

647
	if (!macroFilename.length())
648
	  return;
649

650
651
652
653
654
	// write macro file to the end
	readFileToString(macroFilename, macroData);
	file.seekp(0, ios_base::end);
	long macroPos = file.tellp(), perPos = 0;
	file.write(macroData.c_str(), macroData.length());
655

656
	// write periodic file to the end
657
	if (hasPeriodic) {
658
659
660
661
	  readFileToString(perFilename, periodicData);
	  perPos = file.tellp();
	  file.write(periodicData.c_str(), periodicData.length());
	}
662

663
	// update macroFile_nl
664
665
	int offset = hasPeriodic ? 27 : 13;
	file.seekp(headerLen - offset);
666

667
668
669
670
671
672
	file.write("this:", 5);
	file.write(reinterpret_cast<char*>(&macroPos), 8);
	if (hasPeriodic) {
	  file.write(";this:", 6);
	  file.write(reinterpret_cast<char*>(&perPos), 8);
	}
673
      }
674

675
676
677
678
679
680
681
682
      void readFileToString(std::string filename, std::string& data)
      {
	ifstream file(filename.c_str());
	file.seekg(0, std::ios::end);
	data.reserve(file.tellg());
	file.seekg(0, std::ios::beg);
	data.assign((istreambuf_iterator<char>(file)), istreambuf_iterator<char>());
      }
683
684

      pair<int, int> writeMacroElement(ofstream &file,
Siqi Ling's avatar
Siqi Ling committed
685
686
687
				    MeshStructure &code,
				    vector<vector<double> >& values,
                                    map<const FiniteElemSpace*, vector<int> >& feSpaces,
688
				    Cpsformat::Value cps,
Siqi Ling's avatar
Siqi Ling committed
689
690
691
				    string dataformat)
      {
	stringstream dataStream(ios::out | ios::in | ios::binary);
692

Siqi Ling's avatar
Siqi Ling committed
693
694
        uint32_t nStructureCodes = code.getCode().size();
        dataStream.write(reinterpret_cast<char*>(&nStructureCodes), 4);
695

Siqi Ling's avatar
Siqi Ling committed
696
697
698
        uint32_t codeSize = code.getNumElements();
        dataStream.write(reinterpret_cast<char*>(&codeSize), 4);

699
        dataStream.write(reinterpret_cast<char*>(&(const_cast<vector<uint64_t>&>(code.getCode())[0])),
Siqi Ling's avatar
Siqi Ling committed
700
	       8 * nStructureCodes);
701

Siqi Ling's avatar
Siqi Ling committed
702
703
        if (values.size() > 0) {
	  int moreSize = 0, valuePos = 0;
704
          map<const FiniteElemSpace*, vector<int> >::iterator it;
Siqi Ling's avatar
Siqi Ling committed
705
706
707
708
709
          for(it = feSpaces.begin(); it != feSpaces.end(); it++)
          {
            uint32_t nValuesPerVector = values[valuePos].size();
            dataStream.write(reinterpret_cast<char*>(&nValuesPerVector), 4);
	    moreSize += 4;
710

Siqi Ling's avatar
Siqi Ling committed
711
712
            for (size_t i = 0; i < it->second.size(); i++)
	      moreSize += writeValues(dataStream, dataformat, values[valuePos + i]);
713

Siqi Ling's avatar
Siqi Ling committed
714
715
716
717
718
            valuePos += it->second.size();
          }
        }
	stringstream tmp(ios::out | ios::in | ios::binary);
        boost::iostreams::filtering_streambuf<boost::iostreams::input> in;
719
#ifdef AMDIS_HAS_COMPRESSION
Siqi Ling's avatar
Siqi Ling committed
720
721
	switch(cps)
	{
722
	  case Cpsformat::ZLIB:
Siqi Ling's avatar
Siqi Ling committed
723
724
	    in.push(boost::iostreams::zlib_compressor());
	    break;
725
	  case Cpsformat::BZIP2:
Siqi Ling's avatar
Siqi Ling committed
726
727
	    in.push(boost::iostreams::bzip2_compressor());
	    break;
728
	  case Cpsformat::NONE:
Siqi Ling's avatar
Siqi Ling committed
729
730
731
732
733
734
735
736
737
738
	    break;
	  default:
	    MSG("NOT correct compression flag.\n");
	}
#endif
	in.push(dataStream);
	boost::iostreams::copy(in, tmp);
	file << tmp.rdbuf();
        return make_pair(dataStream.str().length(), tmp.str().length());
      }
739

Siqi Ling's avatar
Siqi Ling committed
740
741
742
      template<typename T>
      int writeValues(stringstream& file, vector<double>& values)
      {
743
744
	size_t size = values.size();
	T* data = new T[size];
Siqi Ling's avatar
Siqi Ling committed
745
746
	for (size_t i = 0; i < size; i++)
	  data[i] = static_cast<T>(values[i]);
747

Siqi Ling's avatar
Siqi Ling committed
748
	file.write(reinterpret_cast<char*>(&data[0]), sizeof(T) * size);
749
	delete [] data;
Siqi Ling's avatar
Siqi Ling committed
750
751
	return sizeof(T) * size;
      }
752

Siqi Ling's avatar
Siqi Ling committed
753
      int writeValues(stringstream& file,
754
		       string dataformat,
Siqi Ling's avatar
Siqi Ling committed
755
756
757
758
759
		       vector<double>& values)
      {
	FUNCNAME("Arh3Writer::detail::writeValues()");
	int newsize = 0;
	std::map<string,Valformat>::const_iterator it = dataformatMap.find(dataformat);
760

Siqi Ling's avatar
Siqi Ling committed
761
	TEST_EXIT(it != dataformatMap.end())("Wrong data format.\n");
762

Siqi Ling's avatar
Siqi Ling committed
763
764
765
766
767
768
769
770
771
	switch(it->second) {
	  case SI08:
	    newsize = writeValues<int8_t>(file, values);
	    break;
	  case SI16:
	    newsize = writeValues<int16_t>(file, values);
	    break;
	  case SI32:
	    newsize = writeValues<int32_t>(file, values);
772
	    break;
Siqi Ling's avatar
Siqi Ling committed
773
774
	  case SI64:
	    newsize = writeValues<int64_t>(file, values);
775
	    break;
Siqi Ling's avatar
Siqi Ling committed
776
777
	  case UI08:
	    newsize = writeValues<uint8_t>(file, values);
778
	    break;
Siqi Ling's avatar
Siqi Ling committed
779
780
781
782
783
	  case UI16:
	    newsize = writeValues<uint16_t>(file, values);
	    break;
	  case UI32:
	    newsize = writeValues<uint32_t>(file, values);
784
	    break;
Siqi Ling's avatar
Siqi Ling committed
785
786
787
788
789
	  case UI64:
	    newsize = writeValues<uint64_t>(file, values);
	    break;
	  case SF32:
	    newsize = writeValues<float>(file, values);
790
	    break;
Siqi Ling's avatar
Siqi Ling committed
791
792
793
794
795
796
797
798
799
800
801
	  case SF64:
	    newsize = writeValues<double>(file, values);
	    break;
	  default:
	    ERROR_EXIT("Wrong data format.\n");
	}
	return newsize;
      }
    }//end namespace detail
  } // end namespace Arh3Writer
} } // end namespace io, AMDiS