ArhReader.cc 12.6 KB
Newer Older
1 2 3
//
// Software License for AMDiS
//
Florian Stenger's avatar
Florian Stenger committed
4
// Copyright (c) 2010 Dresden University of Technology
5 6 7 8 9 10 11 12
// 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
#include <fstream>
14
#include <stdint.h>
15
#include <boost/filesystem.hpp>
16 17 18 19

#include "ArhReader.h"
#include "Mesh.h"
#include "MeshStructure.h"
20 21
#include "Traverse.h"
#include "DOFVector.h"
22
#include "Debug.h"
23 24

namespace AMDiS {
Florian Stenger's avatar
Florian Stenger committed
25

26 27
  using namespace std;

28

29
  void ArhReader::read(string filename, Mesh *mesh,
Florian Stenger's avatar
Florian Stenger committed
30
		       DOFVector<double>* vec0,
31
		       DOFVector<double>* vec1,
32 33 34
		       DOFVector<double>* vec2,
		       bool writeParallel,
		       int nProcs)
35
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
36
    vector<DOFVector<double>*> vecs(0);
37
    if (vec0)
Praetorius, Simon's avatar
Praetorius, Simon committed
38
      vecs.push_back(vec0);
39
    if (vec1)
Praetorius, Simon's avatar
Praetorius, Simon committed
40
      vecs.push_back(vec1);
41
    if (vec2)
Praetorius, Simon's avatar
Praetorius, Simon committed
42
      vecs.push_back(vec2);
Florian Stenger's avatar
Florian Stenger committed
43

44
    ArhReader::read(filename, mesh, vecs, writeParallel, nProcs);
45 46 47
  }


Florian Stenger's avatar
Florian Stenger committed
48
  void ArhReader::read(string filename, Mesh *mesh,
49 50 51
		       vector<DOFVector<double>*> vecs,
		       bool writeParallel,
		       int nProcs)
52 53 54
  {
    FUNCNAME("ArhReader::read()");

55 56 57 58
    if (writeParallel) {
      using boost::lexical_cast;
      int sPos = filename.find(".arh");
      TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n");
Florian Stenger's avatar
Florian Stenger committed
59
      string name = filename.substr(0, sPos);
60

Florian Stenger's avatar
Florian Stenger committed
61
      if (nProcs == -1) {
62
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
63 64
	string procFilename = name + "-p" + lexical_cast<string>(MPI::COMM_WORLD.Get_rank()) + "-.arh";
	readFile(procFilename, mesh, vecs);
65 66 67
#else
	ERROR_EXIT("Reading parallel ARH files in sequential computations requires to specify the number of nodes on which the ARH file was created!\n");
#endif
68 69 70 71 72 73 74 75 76 77 78 79
      } else {
	for (int i = 0; i < nProcs; i++) {
	  string procFilename = name + "-p" + lexical_cast<string>(i) + "-.arh";
	  readFile(procFilename, mesh, vecs);
	}
      }
    } else {
      readFile(filename, mesh, vecs);
    }
    MSG("ARH file read from: %s\n", filename.c_str());
  }

Florian Stenger's avatar
Florian Stenger committed
80

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
  void ArhReader::readMeta(string filename,
			   Mesh *mesh,
			   DOFVector<double>* vec0,
			   DOFVector<double>* vec1,
			   DOFVector<double>* vec2)
  {
    vector<DOFVector<double>*> vecs;
    if (vec0)
      vecs.push_back(vec0);
    if (vec1)
      vecs.push_back(vec1);
    if (vec2)
      vecs.push_back(vec2);

    readMeta(filename, mesh, vecs);
  }


  void ArhReader::readMeta(string filename,
			   Mesh *mesh,
			   vector<DOFVector<double>*> vecs)
  {
    FUNCNAME("ArhReader::readMeta()");

Thomas Witkowski's avatar
Thomas Witkowski committed
105 106 107 108 109 110
    
    // === Read the meta arh file. ===

    string arhPrefix = "";
    map<int, int> elInRank;
    map<int, int> elCodeSize;
111
    int nProc = readMetaData(filename, elInRank, elCodeSize, arhPrefix);
Thomas Witkowski's avatar
Thomas Witkowski committed
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    

    // === Check which arh files must be read by current rank. ===

    // Set of all file indices which should be read to restore all macro elements.
    std::set<int> readArhFiles;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
      int macroElIndex = elInfo->getElement()->getIndex();
      TEST_EXIT(elInRank.count(macroElIndex))("Should not happen!\n");
      readArhFiles.insert(elInRank[macroElIndex]);
      elInfo = stack.traverseNext(elInfo);
    } 


    // === Read the individual arh files. ===

131 132 133
    boost::filesystem::path p(filename.c_str());
    boost::filesystem::path directory = p.parent_path();

Thomas Witkowski's avatar
Thomas Witkowski committed
134 135
    for (std::set<int>::iterator it = readArhFiles.begin();
	 it != readArhFiles.end(); ++it) {
136
      string arhFilename = directory.string() + "/" + arhPrefix;
137 138 139 140 141
      if (nProc == 1)
	arhFilename += ".arh";
      else
	arhFilename += "-p" + boost::lexical_cast<string>(*it) + "-.arh";

Thomas Witkowski's avatar
Thomas Witkowski committed
142 143 144 145 146 147 148 149 150 151 152
      MSG("ARH file read from: %s\n", arhFilename.c_str());
      readFile(arhFilename, mesh, vecs);
    }
  }


  int ArhReader::readMetaData(string filename,
			      map<int, int> &elInRank,
			      map<int, int> &elCodeSize,
			      string &arhPrefix)
  {
153 154 155 156 157 158 159
    ifstream file;
    file.open(filename.c_str());
    TEST_EXIT(file.is_open())
      ("Cannot open arh meta file \"%s\"\n", filename.c_str());

    string readStr = "";
    file >> readStr;
Thomas Witkowski's avatar
Thomas Witkowski committed
160
    arhPrefix = "";
161 162 163 164 165 166 167 168 169 170 171 172
    file >> arhPrefix;
    int nProc;
    file >> nProc;

    // Maps to each macro element index the arh file index it is stored in.
    for (int i = 0; i < nProc; i++) {
      int tmp;
      file >> tmp;
      TEST_EXIT(tmp == i)("Should not happen!\n");
      int nMacroEl;
      file >> nMacroEl;
      for (int j = 0; j < nMacroEl; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
173
	int elIndex, codeSize;
174
	file >> elIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
175 176 177
	file >> codeSize;
	elInRank[elIndex] = i;
	elCodeSize[elIndex] = codeSize;
178 179 180 181 182
      }
    }

    file.close();

Thomas Witkowski's avatar
Thomas Witkowski committed
183 184
    return nProc;
  }
185 186


Thomas Witkowski's avatar
Thomas Witkowski committed
187 188 189
  int ArhReader::readMetaData(string filename)
  {
    FUNCNAME("ArhReader::readMetaData()");
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191 192 193 194 195 196 197 198 199 200 201 202
    ifstream file;
    file.open(filename.c_str());
    TEST_EXIT(file.is_open())
      ("Cannot open arh meta file \"%s\"\n", filename.c_str());

    string readStr = "";
    file >> readStr;
    file >> readStr;
    int nProc;
    file >> nProc;
    file.close();
    return nProc;
203 204 205
  }


Florian Stenger's avatar
Florian Stenger committed
206
  void ArhReader::readFile(string filename,
207 208 209 210 211
			   Mesh *mesh,
			   vector<DOFVector<double>*> vecs)
  {
    FUNCNAME("ArhReader::readFile()");

212 213 214 215 216
    // === Get set of all macro elements in mesh. ===
    std::set<int> macroInMesh;
    for (std::deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
	 it != mesh->getMacroElements().end(); ++it)
      macroInMesh.insert((*it)->getIndex());
Florian Stenger's avatar
Florian Stenger committed
217

218

219
    RefinementManager *refManager = NULL;
220 221 222 223 224 225 226 227 228 229 230
    switch (mesh->getDim()) {
    case 2:
      refManager = new RefinementManager2d();
      break;
    case 3:
      refManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

231

232 233
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
234 235
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
236

237
    string typeId(4, ' ');
238 239 240
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
241 242 243 244 245 246

    file.read(const_cast<char*>(typeId.data()), 4);
    file.read(reinterpret_cast<char*>(&nMacroElements), 4);
    file.read(reinterpret_cast<char*>(&nValueVectors), 4);
    file.read(reinterpret_cast<char*>(&nAllValues), 4);

247
    TEST_EXIT(nValueVectors == vecs.size())
Florian Stenger's avatar
Florian Stenger committed
248
      ("File %s has %d vectors, but %d DOFVectors are provided!\n",
249
       filename.c_str(), nValueVectors, vecs.size());
250

251 252 253 254
    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;
255 256 257 258 259

      file.read(reinterpret_cast<char*>(&elIndex), 4);
      file.read(reinterpret_cast<char*>(&nStructureCodes), 4);
      file.read(reinterpret_cast<char*>(&codeSize), 4);

260
      vector<uint64_t> structureCode(nStructureCodes);
261
      file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);
262 263 264

      MeshStructure elementStructure;
      elementStructure.init(structureCode, codeSize);
Thomas Witkowski's avatar
da  
Thomas Witkowski committed
265
      if (macroInMesh.count(elIndex))
266
	elementStructure.fitMeshToStructure(mesh, refManager, false, elIndex);
267

268
      if (nValueVectors > 0) {
Florian Stenger's avatar
Florian Stenger committed
269 270
	uint32_t nValuesPerVector = 0;
	file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);
271

Florian Stenger's avatar
Florian Stenger committed
272 273 274 275
	for (unsigned int j = 0; j < nValueVectors; j++) {
	  vector<double> values(nValuesPerVector);
	  file.read(reinterpret_cast<char*>(&(values[0])), 8 * nValuesPerVector);
	  if (vecs[j] != NULL) {
276
	    if (macroInMesh.count(elIndex))
Florian Stenger's avatar
Florian Stenger committed
277 278
	      setDofValues(elIndex, mesh, values, vecs[j]);
	  }
279
	}
280 281 282 283 284
      }
    }

    file.close();

285
    delete refManager;
286
  }
287 288 289


  void ArhReader::setDofValues(int macroElIndex, Mesh *mesh,
290
			       vector<double>& values, DOFVector<double>* vec)
291 292 293 294 295 296 297 298 299 300 301
  {
    FUNCNAME("ArhReader::setDofValues()");

    bool macroElement = false;
    int valuePos = 0;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      if (!macroElement) {
	Element *mEl = elInfo->getMacroElement()->getElement();
302
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
Florian Stenger's avatar
Florian Stenger committed
303
	  (*vec)[mEl->getDof(i, 0)] = values[valuePos++];
304 305 306 307 308
	macroElement = true;
      }

      Element *el = elInfo->getElement();
      if (!el->isLeaf())
309
	(*vec)[el->getChild(0)->getDof(mesh->getDim(), 0)] = values[valuePos++];
310 311 312 313

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
314 315


316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
  int ArhReader::getNumValueVectors(string filename)
  {
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
    
    string typeId = "";
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    
    file.read(const_cast<char*>(typeId.data()), 4);
    file.read(reinterpret_cast<char*>(&nMacroElements), 4);
    file.read(reinterpret_cast<char*>(&nValueVectors), 4);
    
    file.close();
    
    return nValueVectors;
  }


Florian Stenger's avatar
Florian Stenger committed
335 336 337
  //the following three functions are identical to read/readfile except that they read from
  //a block of memory instead of from a file
  void ArhReader::readFromMemoryBlock(vector<char> &data, Mesh *mesh,
338 339 340 341 342
				      DOFVector<double>* vec0,
				      DOFVector<double>* vec1,
				      DOFVector<double>* vec2,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
  {
    uint32_t nValueVectors;
    memcpy(reinterpret_cast<char*>(&nValueVectors), &data[8], 4);
    vector<DOFVector<double>*> vecs(0);
    if (nValueVectors > 0)
      vecs.push_back(vec0);
    if (nValueVectors > 1)
      vecs.push_back(vec1);
    if (nValueVectors > 2)
      vecs.push_back(vec2);
    for (uint32_t i = 3; i < nValueVectors; i++)
      vecs.push_back(NULL);

    readFromMemoryBlock(data, mesh, vecs, writeParallel, nProcs);
  }


  void ArhReader::readFromMemoryBlock(vector<char> &data, Mesh *mesh,
361 362 363
				      vector<DOFVector<double>*> vecs,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
364 365 366 367 368 369 370 371 372 373 374 375
  {
    FUNCNAME("ArhReader::readFromMemoryBlock()");

    if (writeParallel) {
      ERROR_EXIT("Reading a hierarchy from several memory blocks is not implemented yet!\n");
    } else {
      readBlock(data, mesh, vecs);
    }
  }


  void ArhReader::readBlock(vector<char> &data,
376 377
			    Mesh *mesh,
			    vector<DOFVector<double>*> vecs)
Florian Stenger's avatar
Florian Stenger committed
378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
  {
    FUNCNAME("ArhReader::readBlock()");

    // === Get set of all macro elements in mesh. ===
    std::set<int> macroInMesh;
    for (std::deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
                                         it != mesh->getMacroElements().end(); ++it)
      macroInMesh.insert((*it)->getIndex());


    RefinementManager *refManager = NULL;
    switch (mesh->getDim()) {
    case 2:
      refManager = new RefinementManager2d();
      break;
    case 3:
      refManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

    string typeId = "";
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
    uint32_t mem_index = 0;

    memcpy(const_cast<char*>(typeId.data()), &data[mem_index], 4);
    mem_index += 4;
    memcpy(reinterpret_cast<char*>(&nMacroElements), &data[mem_index], 4);
    mem_index += 4;
    memcpy(reinterpret_cast<char*>(&nValueVectors), &data[mem_index], 4);
    mem_index += 4;
    memcpy(reinterpret_cast<char*>(&nAllValues), &data[mem_index], 4);
    mem_index += 4;

    TEST_EXIT(nValueVectors == vecs.size())
      ("data has %d vectors, but %d DOFVectors are provided!\n", nValueVectors, vecs.size());

    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;

      memcpy(reinterpret_cast<char*>(&elIndex), &data[mem_index], 4);
      mem_index += 4;
      memcpy(reinterpret_cast<char*>(&nStructureCodes), &data[mem_index], 4);
      mem_index += 4;
      memcpy(reinterpret_cast<char*>(&codeSize), &data[mem_index], 4);
      mem_index += 4;

      vector<uint64_t> structureCode(nStructureCodes);
      memcpy(reinterpret_cast<char*>(&structureCode[0]), &data[mem_index], 8 * nStructureCodes);
      mem_index += 8 * nStructureCodes;

      MeshStructure elementStructure;
      elementStructure.init(structureCode, codeSize);
      if (macroInMesh.count(elIndex) == 1)
        elementStructure.fitMeshToStructure(mesh, refManager, false, elIndex);

      if(nValueVectors > 0){
        uint32_t nValuesPerVector = 0;
        memcpy(reinterpret_cast<char*>(&nValuesPerVector), &data[mem_index], 4);
        mem_index += 4;

        for (unsigned int j = 0; j < nValueVectors; j++) {
          vector<double> values(nValuesPerVector);
          memcpy(reinterpret_cast<char*>(&values[0]), &data[mem_index], 8 * nValuesPerVector);
          mem_index += 8 * nValuesPerVector;
          if (vecs[j] != NULL) {
            if (macroInMesh.count(elIndex) == 1)
              setDofValues(elIndex, mesh, values, vecs[j]);
          }
        }
      }
    }

    delete refManager;
  }

459
}