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
}