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
}