ArhReader.cc 9.55 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 16 17 18

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

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

24 25
  using namespace std;

26

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

42
    ArhReader::read(filename, mesh, vecs, writeParallel, nProcs);
43 44 45
  }


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

53 54 55 56
    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
57
      string name = filename.substr(0, sPos);
58

Florian Stenger's avatar
Florian Stenger committed
59
      if (nProcs == -1) {
60
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
61 62
	string procFilename = name + "-p" + lexical_cast<string>(MPI::COMM_WORLD.Get_rank()) + "-.arh";
	readFile(procFilename, mesh, vecs);
63 64 65
#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
66 67 68 69 70 71 72 73 74 75 76 77
      } 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
78 79

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

85 86 87 88 89
    // === 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
90

91

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

104

105 106
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
107 108
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
109

110
    string typeId(4, ' ');
111 112 113
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
114 115 116 117 118 119

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

120
    TEST_EXIT(nValueVectors == vecs.size())
Florian Stenger's avatar
Florian Stenger committed
121
      ("File %s has %d vectors, but %d DOFVectors are provided!\n",
122
       filename.c_str(), nValueVectors, vecs.size());
123

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

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

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

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

141
      if (nValueVectors > 0) {
Florian Stenger's avatar
Florian Stenger committed
142 143
	uint32_t nValuesPerVector = 0;
	file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);
144

Florian Stenger's avatar
Florian Stenger committed
145 146 147 148 149 150 151
	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) {
	    if (macroInMesh.count(elIndex) == 1)
	      setDofValues(elIndex, mesh, values, vecs[j]);
	  }
152
	}
153 154 155 156 157
      }
    }

    file.close();

158
    delete refManager;
159
  }
160 161 162


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

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

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
187 188


189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
  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
208 209 210
  //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,
211 212 213 214 215
				      DOFVector<double>* vec0,
				      DOFVector<double>* vec1,
				      DOFVector<double>* vec2,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
  {
    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,
234 235 236
				      vector<DOFVector<double>*> vecs,
				      bool writeParallel,
				      int nProcs)
Florian Stenger's avatar
Florian Stenger committed
237 238 239 240 241 242 243 244 245 246 247 248
  {
    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,
249 250
			    Mesh *mesh,
			    vector<DOFVector<double>*> vecs)
Florian Stenger's avatar
Florian Stenger committed
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
  {
    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;
  }

332
}