ProblemImplicit.cc 9.77 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 "ProblemImplicit.h"
Naumann, Andreas's avatar
Naumann, Andreas committed
14
#include "MathFunctions.h"
Naumann, Andreas's avatar
Naumann, Andreas committed
15
#include "io/ArhReader.h"
Naumann, Andreas's avatar
Naumann, Andreas committed
16
using namespace std;
17

18
namespace AMDiS {
19

Naumann, Andreas's avatar
Naumann, Andreas committed
20
  void ProblemImplicit::readDofVec(const std::string& filename, DOFVector<double>* vec, 
21
22
                                       Mesh* mesh) 
  {
Naumann, Andreas's avatar
Naumann, Andreas committed
23
#if 0
24
25
26
27
28
29
    long size = mesh->getNumberOfVertices();

    TEST_EXIT(vec->getSize()>=size)("dofvector smaller than vertex data vector\n");
    std::string buf;
    getline(in, buf);
    while (!in.eof() && buf != "vertex values: solution") 
30
      getline(in, buf);
31
32
33
34
35
36

    TEST_EXIT(!in.eof())("no vertex data\n");

    for (long i = 0; i < size ; ++i) {
      in >> buf;
      (*vec)[i] = atof(buf.c_str());
37
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
38
39
#endif
    ArhReader::read(filename, vec->getFeSpace()->getMesh(), vec);
40
  }
41

42

Naumann, Andreas's avatar
Naumann, Andreas committed
43
  void ProblemImplicit::readR(const std::string& inStream, double eps, 
44
45
	                          Mesh* mesh, int implMesh , int comp) 
  {
Naumann, Andreas's avatar
Naumann, Andreas committed
46
    FUNCNAME("ProblemImplicit::readR()");
47

48
49
50
51
    DOFVector<double>* r = getSignedDistance(implMesh, comp);
    DOFVector<double>* phi1 = getPhi1(implMesh, comp);
    DOFVector<double>* phi2 = getPhi2(implMesh, comp);
    DOFVector<double>* levelSet = getLevelset(implMesh, comp);
52

53
54
55
56
    TEST_EXIT(r != NULL)("no signed distance vector\n");
    TEST_EXIT(phi1 != NULL)("no phasefield1 vector\n");
    TEST_EXIT(phi2 != NULL)("no phasefield2 vector\n");
    TEST_EXIT(levelSet != NULL)("no levelSet vector\n");
57
58

    bool checkSize = r->getSize() == phi1->getSize() && 
59
      r->getSize() == phi2->getSize();
60
61
62
63
64
65
66
67
68

    TEST_EXIT(checkSize)("something went wrong\n");
    
    readDofVec(inStream, r, mesh);

    for (int i = r->getSize() - 1; i >= 0; --i) {
      (*phi1)[i] = Phi1((*r)[i], eps);
      (*phi2)[i] = Phi2((*r)[i], eps);
      (*levelSet)[i] = LevelSet((*r)[i]);
69
    }
70
  }
71

72

Naumann, Andreas's avatar
Naumann, Andreas committed
73
  void ProblemImplicit::readPhi1(const std::string& inStream, double eps, 
74
                                     Mesh* mesh, int implMesh, int comp)
75
  {
76
77
78
79
    DOFVector<double>* r = getSignedDistance(implMesh, comp);
    DOFVector<double>* phi1 = getPhi1(implMesh, comp);
    DOFVector<double>* phi2 = getPhi2(implMesh, comp);
    DOFVector<double>* levelSet = getLevelset(implMesh, comp);
80

81
82
83
84
    TEST_EXIT(r != NULL)("no signed distance vector\n");
    TEST_EXIT(phi1 != NULL)("no phasefield1 vector\n");
    TEST_EXIT(phi2 != NULL)("no phasefield2 vector\n");
    TEST_EXIT(levelSet != NULL)("no levelSet vector\n");
85
86

    bool checkSize = r->getSize() == phi1->getSize() && 
87
      r->getSize() == phi2->getSize();
88
89
90
91
92
93
94
95
96

    TEST_EXIT(checkSize)("something went wrong\n");

    readDofVec(inStream, phi1, mesh);

    for (int i = r->getSize() - 1; i>= 0; --i) {
      (*r)[i] = Phi1ToR((*phi1)[i], eps);
      (*phi2)[i] = 1 - (*phi1)[i];
      (*levelSet)[i] = LevelSet((*r)[i]);
97
    }
98
  }
99

100

Naumann, Andreas's avatar
Naumann, Andreas committed
101
  void ProblemImplicit::readPhi2(const std::string& inStream, double eps, 
102
	                             Mesh* mesh, int implMesh, int comp)
103
  { 
104
105
106
107
    DOFVector<double>* r = getSignedDistance(implMesh, comp);
    DOFVector<double>* phi1 = getPhi1(implMesh, comp);
    DOFVector<double>* phi2 = getPhi2(implMesh, comp);
    DOFVector<double>* levelSet = getLevelset(implMesh, comp);
108

109
110
111
112
    TEST_EXIT(r != NULL)("no signed distance vector\n");
    TEST_EXIT(phi1 != NULL)("no phasefield1 vector\n");
    TEST_EXIT(phi2 != NULL)("no phasefield2 vector\n");
    TEST_EXIT(levelSet != NULL)("no levelSet vector\n");
113
114

    bool checkSize = r->getSize() == phi1->getSize() &&
115
      r->getSize() == phi2->getSize();
116
117
118
119
120
121
122
123
    TEST_EXIT(checkSize)("something went wrong\n");

    readDofVec(inStream, phi2, mesh);

    for (int i = r->getSize() - 1; i >= 0; --i) {
      (*r)[i] = Phi2ToR((*r)[i], eps);
      (*phi1)[i] = 1 - (*phi2)[i];
      (*levelSet)[i] = LevelSet((*r)[i]);
124
125
    }
  }
126
  
127

Naumann, Andreas's avatar
Naumann, Andreas committed
128
  std::string ProblemImplicit::getDofFilename(std::string path, int implMesh)
129
  {
Naumann, Andreas's avatar
Naumann, Andreas committed
130
    FUNCNAME("ProblemImplicit::getDofFilename()");
131
    std::string suffix = "[" + 
132
133
      boost::lexical_cast< std::string >(implMesh) +
      "]";
134
135

    std::string dofFilename("");
136
    Parameters::get(path + "dof file" + suffix, dofFilename);
Naumann, Andreas's avatar
Naumann, Andreas committed
137
138
139
    TEST_EXIT(dofFilename.length() == 0)("dof file was changed to \"arh file\" and reads an arh file now");
    Parameters::get(path + "arh file" + suffix, dofFilename);
    if (implMesh == 0 && dofFilename.length() == 0) {
140
      Parameters::get(path + "dof file", dofFilename);
Naumann, Andreas's avatar
Naumann, Andreas committed
141
142
143
      TEST_EXIT(dofFilename.length() == 0)("dof file was changed to \"arh file\" and reads an arh file now");
      Parameters::get(path + "arh file", dofFilename);
    }
144
145
    return dofFilename;    
  }
146

147

Naumann, Andreas's avatar
Naumann, Andreas committed
148
  double ProblemImplicit::getEpsilon(std::string path, int implMesh)
149
150
  {
    std::string suffix = "[" + 
151
152
      boost::lexical_cast< std::string >(implMesh) +
      "]";
153

154
    double eps(-1);
155
    Parameters::get(path + "eps" + suffix, eps);
156
    if (implMesh == 0 && eps < 0)
157
      Parameters::get(path + "eps", eps);
158
159
160
    return eps;
  }

161

Naumann, Andreas's avatar
Naumann, Andreas committed
162
  int ProblemImplicit::getType(std::string path, int implMesh)
163
  {
164
    std::string suffix = "[" + 
165
166
      boost::lexical_cast< std::string >(implMesh) +
      "]";
167
    int serType(-1);
168
    Parameters::get(path + "type" + suffix, serType);
169
    if (implMesh == 0 && serType < 0)
170
      Parameters::get(path + "type", serType);
171
172
173
174
    return serType;
  }


175
  DOFVector<double>* ProblemImplicit::getSignedDistance(unsigned int im , unsigned int m) 
176
  { 
177
178
179
    if (m >= r.size() || im >= r[m].size())
      return NULL;
    return (r[m])[im]; 
180
181
  }

182

183
  DOFVector<double>* ProblemImplicit::getPhi1(unsigned int im, unsigned int m)
184
  {
185
186
187
188
    if (m >= phi1.size() || im >= phi1[m].size())
      return NULL;

    return (phi1[m])[im];
189
190
  }

191
192

  DOFVector<double>* ProblemImplicit::getPhi2(unsigned int im, unsigned int m)
193
194
195
196
197
198
199
  {
    if (m >= phi2.size() || im >= phi2[m].size())
      return NULL;

    return (phi2[m])[im];
  }

200
201

  DOFVector<double>* ProblemImplicit::getLevelset(unsigned int im, unsigned int m)
202
203
204
205
206
207
  {
    if (m >= levelSet.size() || im >= levelSet[m].size())
      return NULL;

    return (levelSet[m])[im];
  }
208

209
210

  bool ProblemImplicit::createImplicitMesh()
211
212
213
214
215
216
  {
    //check each mesh if it's an implicit one
    r.resize(meshes.size());
    phi1.resize(meshes.size());
    phi2.resize(meshes.size());
    levelSet.resize(meshes.size());
217
    for (unsigned int i = 0; i < meshes.size(); ++i)
218
219
      createImplicitMesh(i);
    return true;
220
221
  }

222

223
  bool ProblemImplicit::createImplicitMesh(int p) 
224
  {
225
226
    FUNCNAME("ProblemImplicit::createImplicitMesh()");

227
    std::string path = name + "->implicit mesh[" 
228
      + boost::lexical_cast< std::string >(p) + "]->";
229
    int nImplMeshes(0);
230
    Parameters::get(path + "nr meshes", nImplMeshes);
Naumann, Andreas's avatar
Naumann, Andreas committed
231
232
233
234
235
236
    if (nImplMeshes == 0)
      return false;
    r[p].resize(nImplMeshes, NULL);
    phi1[p].resize(nImplMeshes, NULL);
    phi2[p].resize(nImplMeshes, NULL);
    levelSet[p].resize(nImplMeshes, NULL);
237
238
239
240
241
242
243
244
245
246
247

    for ( int i = 0; i < nImplMeshes ; ++i ) {
      (r[p])[i] = new DOFVector< double >(getFeSpace(p), "r");
      (phi1[p])[i] = new DOFVector< double >(getFeSpace(p), "phi1");
      (phi2[p])[i] = new DOFVector< double >(getFeSpace(p), "phi2");
      (levelSet[p])[i] = new DOFVector< double >(getFeSpace(p), "levelSet");
      createImplicitMesh(path, i, p);
    }
    return true;
  }
  
248
249
250

  bool ProblemImplicit::createImplicitMesh(std::string path, int implMesh, 
					   int comp)
251
  {
252
253
    FUNCNAME("ProblemImplicit::createImplicitMesh()");

254
    std::string dofFilename = getDofFilename(path, implMesh);
255
256
257
    if (dofFilename.length() == 0)
      return false;

258
    double eps = getEpsilon(path, implMesh);
259
260
    if (eps < 0.0)
      return false;
261
    int serType = getType(path, implMesh);
262
263
    if (serType < 0)
      return false;
264
265
266

    TEST_EXIT(meshes[comp] != NULL)("the mesh was not created\n");
    
267
268
    switch (serType) {
    case 0:
Naumann, Andreas's avatar
Naumann, Andreas committed
269
      readR(dofFilename, eps, meshes[comp], implMesh, comp);
270
271
      break;
    case 1:
Naumann, Andreas's avatar
Naumann, Andreas committed
272
      readPhi1(dofFilename, eps, meshes[comp], implMesh, comp);
273
274
      break;
    case 2:
Naumann, Andreas's avatar
Naumann, Andreas committed
275
      readPhi2(dofFilename, eps, meshes[comp], implMesh, comp);
276
277
278
      break;
    default:
      WARNING("unknown implicit mesh type\n");
279
280
281
282
    }
    return true;
  }

283
284

  void ProblemImplicit::createMesh()
285
  {
286
287
    FUNCNAME("ProblemImplicit::createMesh()");

Naumann, Andreas's avatar
Naumann, Andreas committed
288
289
    ProblemStat::createMesh();
#if 0
290
    implMesh.resize(meshes.size());
291
    for (unsigned int i = 0 ; i < meshes.size() ; ++i ) {
292
      std::string path = name + "->implicit mesh[" +
293
	boost::lexical_cast< std::string >(i) + "]";
294
      std::string meshFilename("");
295
      Parameters::get(path + "->macro file name", meshFilename);
296
      std::string serFilename("");
297
      Parameters::get(path + "->serialization file name", serFilename);
298
299
      implMesh[i] = true;
      if (meshFilename.length() > 0)
300
301
	meshes[i]->setName(path);
      else if (serFilename.length() > 0) {
302
303
304
305
	std::ifstream inStream(serFilename.c_str());
	meshes[i]->deserialize(inStream);
	inStream.close();
      } else
306
	implMesh[i] = false;
Naumann, Andreas's avatar
Naumann, Andreas committed
307
308
    }
#endif    
309
  }
310

311
312
313
314

  void ProblemImplicit::initialize(Flag initFlag, 
				   ProblemStatSeq* adaptProblem, 
				   Flag adoptFlag)
315
  {
316
317
    FUNCNAME("ProblemImplicit::initialize()");

Naumann, Andreas's avatar
Naumann, Andreas committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
    ProblemStat::initialize(initFlag);
    if ( initFlag.isSet(INIT_IMPLICIT_MESH) )
      createImplicitMesh();
  }

  ProblemImplicit::~ProblemImplicit()
  {
    for ( unsigned int p(0); p < meshes.size(); ++p ) {
      for ( unsigned int i = 0; i < r[p].size() ; ++i )
	  if ( r[p][i] != NULL)
	    delete r[p][i];
       for ( unsigned int i(0); i < phi1[p].size(); ++i )
	 if ( phi1[p][i] != NULL)
	    delete phi1[p][i];
       for ( unsigned int i(0); i < phi2[p].size(); ++i )
	 if ( phi2[p][i] != NULL)
	    delete phi2[p][i];
       for ( unsigned int i(0); i < levelSet[p].size(); ++i )
	 if ( levelSet[p][i] != NULL)
	    delete levelSet[p][i];
    }

340
341
  }
}