ArhReader.cc 12.5 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 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    
    // === Read the meta arh file. ===

    string arhPrefix = "";
    map<int, int> elInRank;
    map<int, int> elCodeSize;
    readMetaData(filename, elInRank, elCodeSize, arhPrefix);
    

    // === 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 136 137 138 139 140 141 142 143 144 145 146 147 148
    for (std::set<int>::iterator it = readArhFiles.begin();
	 it != readArhFiles.end(); ++it) {
      string arhFilename = 
	directory.native() + "/" + arhPrefix + "-p" + boost::lexical_cast<string>(*it) + "-.arh";
      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)
  {
149 150 151 152 153 154 155
    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
156
    arhPrefix = "";
157 158 159 160 161 162 163 164 165 166 167 168
    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
169
	int elIndex, codeSize;
170
	file >> elIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
171 172 173
	file >> codeSize;
	elInRank[elIndex] = i;
	elCodeSize[elIndex] = codeSize;
174 175 176 177 178
      }
    }

    file.close();

Thomas Witkowski's avatar
Thomas Witkowski committed
179 180
    return nProc;
  }
181 182


Thomas Witkowski's avatar
Thomas Witkowski committed
183 184 185
  int ArhReader::readMetaData(string filename)
  {
    FUNCNAME("ArhReader::readMetaData()");
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187 188 189 190 191 192 193 194 195 196 197 198
    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;
199 200 201
  }


Florian Stenger's avatar
Florian Stenger committed
202
  void ArhReader::readFile(string filename,
203 204 205 206 207
			   Mesh *mesh,
			   vector<DOFVector<double>*> vecs)
  {
    FUNCNAME("ArhReader::readFile()");

208 209 210 211 212
    // === 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
213

214

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

227

228 229
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
230 231
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
232

233
    string typeId(4, ' ');
234 235 236
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
237 238 239 240 241 242

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

243
    TEST_EXIT(nValueVectors == vecs.size())
Florian Stenger's avatar
Florian Stenger committed
244
      ("File %s has %d vectors, but %d DOFVectors are provided!\n",
245
       filename.c_str(), nValueVectors, vecs.size());
246

247 248 249 250
    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;
251 252 253 254 255

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

256
      vector<uint64_t> structureCode(nStructureCodes);
257
      file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);
258 259 260

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

264
      if (nValueVectors > 0) {
Florian Stenger's avatar
Florian Stenger committed
265 266
	uint32_t nValuesPerVector = 0;
	file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);
267

Florian Stenger's avatar
Florian Stenger committed
268 269 270 271
	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) {
272
	    if (macroInMesh.count(elIndex))
Florian Stenger's avatar
Florian Stenger committed
273 274
	      setDofValues(elIndex, mesh, values, vecs[j]);
	  }
275
	}
276 277 278 279 280
      }
    }

    file.close();

281
    delete refManager;
282
  }
283 284 285


  void ArhReader::setDofValues(int macroElIndex, Mesh *mesh,
286
			       vector<double>& values, DOFVector<double>* vec)
287 288 289 290 291 292 293 294 295 296 297
  {
    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();
298
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
Florian Stenger's avatar
Florian Stenger committed
299
	  (*vec)[mEl->getDof(i, 0)] = values[valuePos++];
300 301 302 303 304
	macroElement = true;
      }

      Element *el = elInfo->getElement();
      if (!el->isLeaf())
305
	(*vec)[el->getChild(0)->getDof(mesh->getDim(), 0)] = values[valuePos++];
306 307 308 309

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
310 311


312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
  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
331 332 333
  //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,
334 335 336 337 338
				      DOFVector<double>* vec0,
				      DOFVector<double>* vec1,
				      DOFVector<double>* vec2,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
  {
    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,
357 358 359
				      vector<DOFVector<double>*> vecs,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
360 361 362 363 364 365 366 367 368 369 370 371
  {
    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,
372 373
			    Mesh *mesh,
			    vector<DOFVector<double>*> vecs)
Florian Stenger's avatar
Florian Stenger committed
374 375 376 377 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
  {
    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;
  }

455
}