ArhReader.cc 9.28 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
    size_t 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);
Florian Stenger's avatar
Florian Stenger committed
44

45
    ArhReader::read(filename, mesh, vecs, writeParallel, nProcs);
46 47 48
  }


Florian Stenger's avatar
Florian Stenger committed
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
    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
60
      string name = filename.substr(0, sPos);
61

Florian Stenger's avatar
Florian Stenger committed
62
      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
      } 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
81 82

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

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

94

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
    ifstream file;
    file.open(filename.c_str(), ios::in | ios::binary);
110 111
    TEST_EXIT(file.is_open())
      ("Cannot open file %s\n", filename.c_str());
112

113
    string typeId(4, ' ');
114 115 116
    uint32_t nMacroElements = 0;
    uint32_t nValueVectors = 0;
    uint32_t nAllValues = 0;
117 118 119 120 121 122

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

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

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

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

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

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

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

Florian Stenger's avatar
Florian Stenger committed
148 149 150 151 152 153 154
	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]);
	  }
155
	}
156 157 158 159 160
      }
    }

    file.close();

161
    delete refManager;
162
  }
163 164 165


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

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

      elInfo = stack.traverseNext(elInfo);
    }
  }
Florian Stenger's avatar
Florian Stenger committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 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


  //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,
           DOFVector<double>* vec0,
           DOFVector<double>* vec1,
           DOFVector<double>* vec2,
           bool writeParallel,
           int nProcs)
  {
    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,
           vector<DOFVector<double>*> vecs,
           bool writeParallel,
           int nProcs)
  {
    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,
         Mesh *mesh,
         vector<DOFVector<double>*> vecs)
  {
    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;
  }

316
}