ArhReader.cc 13 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
#include <fstream>
23
#include <stdint.h>
24
#include <boost/filesystem.hpp>
25 26 27 28

#include "ArhReader.h"
#include "Mesh.h"
#include "MeshStructure.h"
29 30
#include "Traverse.h"
#include "DOFVector.h"
31
#include "Debug.h"
32 33

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

35 36
  using namespace std;

37

38
  void ArhReader::read(string filename, Mesh *mesh,
Florian Stenger's avatar
Florian Stenger committed
39
		       DOFVector<double>* vec0,
40
		       DOFVector<double>* vec1,
41 42 43
		       DOFVector<double>* vec2,
		       bool writeParallel,
		       int nProcs)
44
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
45
    vector<DOFVector<double>*> vecs(0);
46
    if (vec0)
Praetorius, Simon's avatar
Praetorius, Simon committed
47
      vecs.push_back(vec0);
48
    if (vec1)
Praetorius, Simon's avatar
Praetorius, Simon committed
49
      vecs.push_back(vec1);
50
    if (vec2)
Praetorius, Simon's avatar
Praetorius, Simon committed
51
      vecs.push_back(vec2);
Florian Stenger's avatar
Florian Stenger committed
52

53
    ArhReader::read(filename, mesh, vecs, writeParallel, nProcs);
54 55 56
  }


Florian Stenger's avatar
Florian Stenger committed
57
  void ArhReader::read(string filename, Mesh *mesh,
58 59 60
		       vector<DOFVector<double>*> vecs,
		       bool writeParallel,
		       int nProcs)
61 62 63
  {
    FUNCNAME("ArhReader::read()");

64 65 66 67
    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
68
      string name = filename.substr(0, sPos);
69

Florian Stenger's avatar
Florian Stenger committed
70
      if (nProcs == -1) {
71
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
72 73
	string procFilename = name + "-p" + lexical_cast<string>(MPI::COMM_WORLD.Get_rank()) + "-.arh";
	readFile(procFilename, mesh, vecs);
74 75 76
#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
77 78 79 80 81 82 83 84 85 86 87 88
      } 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
89

90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
  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
114 115 116 117 118 119
    
    // === Read the meta arh file. ===

    string arhPrefix = "";
    map<int, int> elInRank;
    map<int, int> elCodeSize;
120
    int nProc = readMetaData(filename, elInRank, elCodeSize, arhPrefix);
Thomas Witkowski's avatar
Thomas Witkowski committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
    

    // === 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. ===

140 141 142
    boost::filesystem::path p(filename.c_str());
    boost::filesystem::path directory = p.parent_path();

Thomas Witkowski's avatar
Thomas Witkowski committed
143 144
    for (std::set<int>::iterator it = readArhFiles.begin();
	 it != readArhFiles.end(); ++it) {
145
      string arhFilename = directory.string() + "/" + arhPrefix;
146 147 148 149 150
      if (nProc == 1)
	arhFilename += ".arh";
      else
	arhFilename += "-p" + boost::lexical_cast<string>(*it) + "-.arh";

Thomas Witkowski's avatar
Thomas Witkowski committed
151 152 153 154 155 156 157 158 159 160 161
      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)
  {
162 163 164 165 166 167 168
    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
169
    arhPrefix = "";
170 171 172 173 174 175 176 177 178 179 180 181
    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
182
	int elIndex, codeSize;
183
	file >> elIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
184 185 186
	file >> codeSize;
	elInRank[elIndex] = i;
	elCodeSize[elIndex] = codeSize;
187 188 189 190 191
      }
    }

    file.close();

Thomas Witkowski's avatar
Thomas Witkowski committed
192 193
    return nProc;
  }
194 195


Thomas Witkowski's avatar
Thomas Witkowski committed
196 197 198
  int ArhReader::readMetaData(string filename)
  {
    FUNCNAME("ArhReader::readMetaData()");
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200 201 202 203 204 205 206 207 208 209 210 211
    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;
212 213 214
  }


Florian Stenger's avatar
Florian Stenger committed
215
  void ArhReader::readFile(string filename,
216 217 218 219 220
			   Mesh *mesh,
			   vector<DOFVector<double>*> vecs)
  {
    FUNCNAME("ArhReader::readFile()");

221 222 223 224 225
    // === 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
226

227

228
    RefinementManager *refManager = nullptr;
229 230 231 232 233 234 235 236 237 238 239
    switch (mesh->getDim()) {
    case 2:
      refManager = new RefinementManager2d();
      break;
    case 3:
      refManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

240

241 242
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
243 244
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
245

246
    string typeId(4, ' ');
247 248 249
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
250 251 252 253 254 255

    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);

256
    TEST_EXIT(nValueVectors == vecs.size())
Florian Stenger's avatar
Florian Stenger committed
257
      ("File %s has %d vectors, but %d DOFVectors are provided!\n",
258
       filename.c_str(), nValueVectors, vecs.size());
259

260 261 262 263
    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;
264 265 266 267 268

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

269
      vector<uint64_t> structureCode(nStructureCodes);
270
      file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);
271 272 273

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

277
      if (nValueVectors > 0) {
Florian Stenger's avatar
Florian Stenger committed
278 279
	uint32_t nValuesPerVector = 0;
	file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);
280

Florian Stenger's avatar
Florian Stenger committed
281 282 283
	for (unsigned int j = 0; j < nValueVectors; j++) {
	  vector<double> values(nValuesPerVector);
	  file.read(reinterpret_cast<char*>(&(values[0])), 8 * nValuesPerVector);
284
	  if (vecs[j] != nullptr) {
285
	    if (macroInMesh.count(elIndex))
Florian Stenger's avatar
Florian Stenger committed
286 287
	      setDofValues(elIndex, mesh, values, vecs[j]);
	  }
288
	}
289 290 291 292 293
      }
    }

    file.close();

294
    delete refManager;
295
  }
296 297 298


  void ArhReader::setDofValues(int macroElIndex, Mesh *mesh,
299
			       vector<double>& values, DOFVector<double>* vec)
300 301 302 303 304 305 306 307 308
  {
    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();
309
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
Florian Stenger's avatar
Florian Stenger committed
310
	  (*vec)[mEl->getDof(i, 0)] = values[valuePos++];
311 312 313 314 315
	macroElement = true;
      }

      Element *el = elInfo->getElement();
      if (!el->isLeaf())
316
	(*vec)[el->getChild(0)->getDof(mesh->getDim(), 0)] = values[valuePos++];
317 318 319 320

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
321 322


323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
  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
342 343 344
  //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,
345 346 347 348 349
				      DOFVector<double>* vec0,
				      DOFVector<double>* vec1,
				      DOFVector<double>* vec2,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
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++)
361
      vecs.push_back(nullptr);
Florian Stenger's avatar
Florian Stenger committed
362 363 364 365 366 367

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


  void ArhReader::readFromMemoryBlock(vector<char> &data, Mesh *mesh,
368 369 370
				      vector<DOFVector<double>*> vecs,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
371 372 373 374 375 376 377 378 379 380 381 382
  {
    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,
383 384
			    Mesh *mesh,
			    vector<DOFVector<double>*> vecs)
Florian Stenger's avatar
Florian Stenger committed
385 386 387 388 389 390 391 392 393 394
  {
    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());


395
    RefinementManager *refManager = nullptr;
Florian Stenger's avatar
Florian Stenger committed
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
    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;
455
          if (vecs[j] != nullptr) {
Florian Stenger's avatar
Florian Stenger committed
456 457 458 459 460 461 462 463 464 465
            if (macroInMesh.count(elIndex) == 1)
              setDofValues(elIndex, mesh, values, vecs[j]);
          }
        }
      }
    }

    delete refManager;
  }

466
}