ProblemImplicit.cc 12.1 KB
Newer Older
1
#include "ProblemImplicit.h"
Naumann, Andreas's avatar
Naumann, Andreas committed
2
3
4
#include "MathFunctions.h"

using namespace std;
5
namespace AMDiS {
6
7
8
9
10
11
12
13
14
  void ProblemImplicitBase::readDofVec(std::istream& in, DOFVector<double>* vec, 
                                       Mesh* mesh) 
  {
    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") 
15
      getline(in, buf);
16
17
18
19
20
21

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

    for (long i = 0; i < size ; ++i) {
      in >> buf;
      (*vec)[i] = atof(buf.c_str());
22
    }
23
  }
24

25
26
27
28
  void ProblemImplicitBase::readR(std::istream& inStream, double eps, 
	                          Mesh* mesh, int implMesh , int comp) 
  {
    FUNCNAME("ProblemImplicitBase::readR()");
29

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
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);

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

    bool checkSize = r->getSize() == phi1->getSize() && 
	             r->getSize() == phi2->getSize();

    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]);
52
    }
53
  }
54

55

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
  void ProblemImplicitBase::readPhi1(istream& inStream, double eps, 
                                     Mesh* mesh, int implMesh, int comp )
  {
    DOFVector< double >* r = getSignedDistance(implMesh, comp);
    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
    DOFVector< double >* levelSet = getLevelset(implMesh, comp);

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

    bool checkSize = r->getSize() == phi1->getSize() && 
	             r->getSize() == phi2->getSize();

    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]);
80
    }
81
  }
82

83

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
  void ProblemImplicitBase::readPhi2(istream& inStream, double eps, 
	                             Mesh* mesh, int implMesh, int comp )
  { 
    DOFVector< double >* r = getSignedDistance(implMesh, comp);
    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
    DOFVector< double >* levelSet = getLevelset(implMesh, comp);

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

    bool checkSize = r->getSize() == phi1->getSize() &&
	             r->getSize() == phi2->getSize();
    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]);
107
108
    }
  }
109
110
111
112
113
114
115
116
117
118
119
120
121
  
  std::string ProblemImplicitBase::getDofFilename(std::string path, int implMesh)
  {
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";

    std::string dofFilename("");
    GET_PARAMETER(0, path + "dof file" + suffix, &dofFilename);
    if (implMesh == 0 && dofFilename.length() == 0)
      GET_PARAMETER(0, path + "dof file", &dofFilename);
    return dofFilename;    
  }
122

123
124
125
126
127
  double ProblemImplicitBase::getEpsilon(std::string path, int implMesh)
  {
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";
128

129
130
131
132
133
134
135
136
    double eps(-1);
    GET_PARAMETER(0, path + "eps" + suffix, "%g", &eps);
    if (implMesh == 0 && eps < 0)
      GET_PARAMETER(0, path + "eps", "%g", &eps);
    return eps;
  }

  int ProblemImplicitBase::getType(std::string path, int implMesh)
137
  {
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";
    int serType(-1);
    GET_PARAMETER(0, path + "type" + suffix, "%d", &serType);
    if (implMesh == 0 && serType < 0)
      GET_PARAMETER(0, path + "type", "%d", &serType);
    return serType;
  }

  DOFVector< double >* ProblemImplicitScal::getSignedDistance(int im , int )
  { 
    if (readImplMesh)
      return r[im]; 
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getPhi1(int im , int )
  { 
    if (readImplMesh)
      return phi1[im];
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getPhi2(int im , int )
  { 
    if (readImplMesh)
      return phi2[im]; 
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getLevelset(int im , int ) 
  { 
    if (readImplMesh)
      return levelSet[im]; 
    return NULL;
  }

  void ProblemImplicitScal::createMesh()
  {
    FUNCNAME("ProblemImplicitScal::createMesh()");
    ProblemScal::createMesh();
    std::string path = name + "->implicit mesh";
181
    std::string meshFilename("");
182
183
184
185

    GET_PARAMETER(0, path + "->macro file name", &meshFilename);
    std::string serFilename("");
    GET_PARAMETER(0, path + "->serialization file name", &serFilename);
186
    if ( meshFilename.length() > 0) {
187
      mesh->setName(path);
188
189
      readImplMesh = true;
    }
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    else 
      if ( serFilename.length() > 0 ) {
        //INFO(5,1)("loading implicit mesh %s \n", serFilename.c_str());
	std::string oldName = mesh->getName();
        std::ifstream inStream(serFilename.c_str());
        mesh->deserialize(inStream);
        inStream.close();
	mesh->setName(oldName);
        readImplMesh = true;
      }

  }

  void ProblemImplicitScal::initialize(Flag initFlag, ProblemScal* adaptProblem, 
		                       Flag adoptFlag)
  {
     FUNCNAME("ProblemImplicitScal::initialize()");
     ProblemScal::initialize(initFlag);
     createImplicitMesh();
  }

  bool ProblemImplicitScal::createImplicitMesh() 
  {
213
214
    FUNCNAME("ProblemImplicitScal::createImplicitMesh()");
    INFO(5,1)("setting implicit mesh\n");
215
    if (!isImplicitMesh())
216
      return false;
217
218
219
    std::string path = name + "->implicit mesh->";
    int nImplMesh(1);
    GET_PARAMETER(0, name + "nr meshes", "%d", &nImplMesh);
220

221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
    r.resize(nImplMesh);
    phi1.resize(nImplMesh);
    phi2.resize(nImplMesh);
    levelSet.resize(nImplMesh);
    for ( int i = 0; i < nImplMesh ; ++i) {
      r[i] = new DOFVector< double >(getFeSpace(), "r");
      phi1[i] = new DOFVector< double >(getFeSpace(), "phi1");
      phi2[i] = new DOFVector< double >(getFeSpace(), "phi2");
      levelSet[i] = new DOFVector< double >(getFeSpace(), "levelSet");
      createImplicitMesh(path, i);
    }
    return true;
  }

  bool ProblemImplicitScal::createImplicitMesh(string path, int p)
  {
    
    std::string dofFilename = getDofFilename(path, p);
239
240
241
    if (dofFilename.length() == 0)
      return false;

242
243
    double eps = getEpsilon(path, p);
    if (eps < 0.0)
244
      return false;
245
    int serType = getType(path, p);
246
    if (serType < 0)
247
248
249
      return false;

    TEST_EXIT(mesh != NULL)("the mesh was not created\n");
250
251
    
    std::ifstream inStream(dofFilename.c_str());
252
253
    switch (serType) {
    case 0:
254
      readR(inStream, eps, mesh, p);
255
256
      break;
    case 1:
257
      readPhi1(inStream, eps, mesh, p);
258
259
      break;
    case 2:
260
      readPhi2(inStream, eps, mesh, p);
261
262
263
264
265
266
267
268
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

269

270
  DOFVector< double >* ProblemImplicitVec::getSignedDistance(unsigned int im , unsigned int m) 
271
272
273
274
275
276
  { 
    if ( m >= r.size() || im >= r[m].size())
	return NULL;
      return (r[m])[im]; 
    }

277
  DOFVector< double >* ProblemImplicitVec::getPhi1(unsigned int im, unsigned int m)
278
  {
279
280
281
282
    if (m >= phi1.size() || im >= phi1[m].size())
      return NULL;

    return (phi1[m])[im];
283
284
  }

285
  DOFVector< double >* ProblemImplicitVec::getPhi2(unsigned int im, unsigned int m)
286
287
288
289
290
291
292
  {
    if (m >= phi2.size() || im >= phi2[m].size())
      return NULL;

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

293
  DOFVector< double >* ProblemImplicitVec::getLevelset(unsigned int im, unsigned int m)
294
295
296
297
298
299
  {
    if (m >= levelSet.size() || im >= levelSet[m].size())
      return NULL;

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

301
302
303
304
305
306
307
  bool ProblemImplicitVec::createImplicitMesh()
  {
    //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());
308
    for (unsigned int i = 0; i < meshes.size(); ++i)
309
310
      createImplicitMesh(i);
    return true;
311
312
  }

313

314
315
  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
316
317
318
    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
    //INFO(5,1)("initing mesh %d\n", p);
    if (!isImplicitMesh(p))
319
      return false;
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
    std::string path = name + "->implicit mesh[" 
	             + boost::lexical_cast< std::string >(p) + "]->";
    int nImplMeshes(0);
    GET_PARAMETER(0, path + "nr meshes", "%d", &nImplMeshes);
    //INFO(5,1)("has %d meshe(s)\n", nImplMeshes);
    r[p].resize(nImplMeshes);
    phi1[p].resize(nImplMeshes);
    phi2[p].resize(nImplMeshes);
    levelSet[p].resize(nImplMeshes);

    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;
  }
  
  bool ProblemImplicitVec::createImplicitMesh(std::string path, int implMesh, 
		                              int comp)
  {
    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
    std::string dofFilename = getDofFilename(path, implMesh);
345
346
347
    if (dofFilename.length() == 0)
      return false;

348
    double eps = getEpsilon(path, implMesh);
349
350
    if (eps < 0.0)
      return false;
351
    int serType = getType(path, implMesh);
352
353
    if (serType < 0)
      return false;
354
355
356
357

    TEST_EXIT(meshes[comp] != NULL)("the mesh was not created\n");
    
    std::ifstream inStream(dofFilename.c_str());
358
359
    switch (serType) {
    case 0:
360
      readR(inStream, eps, meshes[comp], implMesh, comp);
361
362
      break;
    case 1:
363
      readPhi1(inStream, eps, meshes[comp], implMesh, comp);
364
365
      break;
    case 2:
366
      readPhi2(inStream, eps, meshes[comp], implMesh, comp);
367
368
369
      break;
    default:
      WARNING("unknown implicit mesh type\n");
370
371
372
373
374
    }
    inStream.close();
    return true;
  }

375
376
377
378
379
  void ProblemImplicitVec::createMesh()
  {
    FUNCNAME("ProblemImplicitVec::createMesh()");
    ProblemVec::createMesh();
    implMesh.resize(meshes.size());
380
    for (unsigned int i = 0 ; i < meshes.size() ; ++i )
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
    {
      std::string path = name + "->implicit mesh[" +
                         boost::lexical_cast< std::string >(i) + "]";
      std::string meshFilename("");
      GET_PARAMETER(0, path + "->macro file name", &meshFilename);
      std::string serFilename("");
      GET_PARAMETER(0, path + "->serialization file name", &serFilename);
      implMesh[i] = true;
      if (meshFilename.length() > 0)
        meshes[i]->setName(path);
      else if (serFilename.length() > 0) 
      {
        //INFO(5,1)("reading serialization file name %s \n", serFilename.c_str());
	std::ifstream inStream(serFilename.c_str());
	meshes[i]->deserialize(inStream);
	inStream.close();
      } else
        implMesh[i] = false;
    }

  }
402

403
404
  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, 
		                      Flag adoptFlag)
405
  {
406
     FUNCNAME("ProblemImplicitVec::initialize()");
407
     ProblemVec::initialize(initFlag);    
408
     createImplicitMesh();
409
410
  }
}