ArhReader.cc 11.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 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 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
  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()");

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

    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;
    string arhPrefix = "";
    file >> arhPrefix;
    int nProc;
    file >> nProc;

    // Maps to each macro element index the arh file index it is stored in.
    map<int, int> macroInProc;
    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++) {
	int elIndex;
	file >> elIndex;
	macroInProc[elIndex] = i;
      }
    }

    file.close();


    // 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(macroInProc.count(macroElIndex))("Should not happen!\n");
      readArhFiles.insert(macroInProc[macroElIndex]);
      elInfo = stack.traverseNext(elInfo);
    } 

    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("READ FILE: %s\n", arhFilename.c_str());
      readFile(arhFilename, mesh, vecs);
    }
  }


Florian Stenger's avatar
Florian Stenger committed
160
  void ArhReader::readFile(string filename,
161 162 163 164 165
			   Mesh *mesh,
			   vector<DOFVector<double>*> vecs)
  {
    FUNCNAME("ArhReader::readFile()");

166 167 168 169 170
    // === 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
171

172

173
    RefinementManager *refManager = NULL;
174 175 176 177 178 179 180 181 182 183 184
    switch (mesh->getDim()) {
    case 2:
      refManager = new RefinementManager2d();
      break;
    case 3:
      refManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

185

186 187
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
188 189
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
190

191
    string typeId(4, ' ');
192 193 194
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
195 196 197 198 199 200

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

201
    TEST_EXIT(nValueVectors == vecs.size())
Florian Stenger's avatar
Florian Stenger committed
202
      ("File %s has %d vectors, but %d DOFVectors are provided!\n",
203
       filename.c_str(), nValueVectors, vecs.size());
204

205 206 207 208
    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;
209 210 211 212 213

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

214
      vector<uint64_t> structureCode(nStructureCodes);
215
      file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);
216 217 218

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

222
      if (nValueVectors > 0) {
Florian Stenger's avatar
Florian Stenger committed
223 224
	uint32_t nValuesPerVector = 0;
	file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);
225

Florian Stenger's avatar
Florian Stenger committed
226 227 228 229
	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) {
230
	    if (macroInMesh.count(elIndex))
Florian Stenger's avatar
Florian Stenger committed
231 232
	      setDofValues(elIndex, mesh, values, vecs[j]);
	  }
233
	}
234 235 236 237 238
      }
    }

    file.close();

239
    delete refManager;
240
  }
241 242 243


  void ArhReader::setDofValues(int macroElIndex, Mesh *mesh,
244
			       vector<double>& values, DOFVector<double>* vec)
245 246 247 248 249 250 251 252 253 254 255
  {
    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();
256
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
Florian Stenger's avatar
Florian Stenger committed
257
	  (*vec)[mEl->getDof(i, 0)] = values[valuePos++];
258 259 260 261 262
	macroElement = true;
      }

      Element *el = elInfo->getElement();
      if (!el->isLeaf())
263
	(*vec)[el->getChild(0)->getDof(mesh->getDim(), 0)] = values[valuePos++];
264 265 266 267

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
268 269


270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
  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
289 290 291
  //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,
292 293 294 295 296
				      DOFVector<double>* vec0,
				      DOFVector<double>* vec1,
				      DOFVector<double>* vec2,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
  {
    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,
315 316 317
				      vector<DOFVector<double>*> vecs,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
318 319 320 321 322 323 324 325 326 327 328 329
  {
    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,
330 331
			    Mesh *mesh,
			    vector<DOFVector<double>*> vecs)
Florian Stenger's avatar
Florian Stenger committed
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 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
  {
    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;
  }

413
}