ArhReader.cc 5.18 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// 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 16 17 18

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

namespace AMDiS {
  
  using namespace std;

26

27 28 29
  void ArhReader::read(string filename, Mesh *mesh,
		       DOFVector<double>* vec0, 
		       DOFVector<double>* vec1,
30 31 32
		       DOFVector<double>* vec2,
		       bool writeParallel,
		       int nProcs)
33
  {
34
    int nValueVectors = getNumValueVectors(filename);
Praetorius, Simon's avatar
Praetorius, Simon committed
35
    vector<DOFVector<double>*> vecs(0);
36
    if (nValueVectors > 0)
Praetorius, Simon's avatar
Praetorius, Simon committed
37
      vecs.push_back(vec0);
38
    if (nValueVectors > 1)
Praetorius, Simon's avatar
Praetorius, Simon committed
39
      vecs.push_back(vec1);
40
    if (nValueVectors > 2)
Praetorius, Simon's avatar
Praetorius, Simon committed
41 42 43
      vecs.push_back(vec2);
    for (size_t i = 3; i < nValueVectors; i++)
      vecs.push_back(NULL);
44
    
45
    ArhReader::read(filename, mesh, vecs, writeParallel, nProcs);
46 47 48 49
  }


  void ArhReader::read(string filename, Mesh *mesh, 
50 51 52
		       vector<DOFVector<double>*> vecs,
		       bool writeParallel,
		       int nProcs)
53 54 55
  {
    FUNCNAME("ArhReader::read()");

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

      if (nProcs == -1) {     
63
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
64 65
	string procFilename = name + "-p" + lexical_cast<string>(MPI::COMM_WORLD.Get_rank()) + "-.arh";
	readFile(procFilename, mesh, vecs);
66 67 68
#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
69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
      } 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());
  }

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

88 89 90 91 92 93 94
    // === 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());
    

95
    RefinementManager *refManager = NULL;
96 97 98 99 100 101 102 103 104 105 106
    switch (mesh->getDim()) {
    case 2:
      refManager = new RefinementManager2d();
      break;
    case 3:
      refManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

107

108 109 110 111
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);

    string typeId = "";
112 113 114
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
115 116 117 118 119 120

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

121 122 123 124
    TEST_EXIT(nValueVectors == vecs.size())
      ("File has %d vectors, but %d DOFVectors are provided!\n", 
       nValueVectors, vecs.size());

125 126 127 128
    for (unsigned int i = 0; i < nMacroElements; i++) {
      uint32_t elIndex = 0;
      uint32_t nStructureCodes = 0;
      uint32_t codeSize = 0;
129 130 131 132 133

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

134
      vector<uint64_t> structureCode(nStructureCodes);
135
      file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);
136 137 138

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

142
      uint32_t nValuesPerVector = 0;
143 144
      file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);

145 146 147
      for (unsigned int j = 0; j < nValueVectors; j++) {
	vector<double> values(nValuesPerVector);
	file.read(reinterpret_cast<char*>(&(values[0])), 8 * nValuesPerVector);
148 149 150 151
	if (vecs[j] != NULL) {
	  if (macroInMesh.count(elIndex) == 1)
	    setDofValues(elIndex, mesh, values, vecs[j]);
	}
152 153 154 155 156
      }
    }

    file.close();

157
    delete refManager;
158
  }
159 160 161


  void ArhReader::setDofValues(int macroElIndex, Mesh *mesh,
162
			       vector<double>& values, DOFVector<double>* vec)
163 164 165 166 167 168 169 170 171 172 173
  {
    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();
174
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
175
	  (*vec)[mEl->getDof(i, 0)] = values[valuePos++];	
176 177 178 179 180
	macroElement = true;
      }

      Element *el = elInfo->getElement();
      if (!el->isLeaf())
181
	(*vec)[el->getChild(0)->getDof(mesh->getDim(), 0)] = values[valuePos++];
182 183 184 185

      elInfo = stack.traverseNext(elInfo);
    }
  }
186
}