ProblemImplicit.cc 11.8 KB
Newer Older
1
2
3
#include "ProblemImplicit.h"

namespace AMDiS {
4
5
6
7
8
9
10
11
12
  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") 
13
      getline(in, buf);
14
15
16
17
18
19

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

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

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

28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
    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]);
50
    }
51
  }
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
  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]);
78
    }
79
  }
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
  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]);
105
106
    }
  }
107
108
109
110
111
112
113
114
115
116
117
118
119
  
  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;    
  }
120

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

127
128
129
130
131
132
133
134
    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)
135
  {
136
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
    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";
179
    std::string meshFilename("");
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

    GET_PARAMETER(0, path + "->macro file name", &meshFilename);
    std::string serFilename("");
    GET_PARAMETER(0, path + "->serialization file name", &serFilename);
    if ( meshFilename.length() > 0)
      mesh->setName(path);
    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() 
  {
    if (!isImplicitMesh())
210
      return false;
211
212
213
    std::string path = name + "->implicit mesh->";
    int nImplMesh(1);
    GET_PARAMETER(0, name + "nr meshes", "%d", &nImplMesh);
214

215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    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);
233
234
235
    if (dofFilename.length() == 0)
      return false;

236
237
    double eps = getEpsilon(path, p);
    if (eps < 0.0)
238
      return false;
239
    int serType = getType(path, p);
240
    if (serType < 0)
241
242
243
      return false;

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

263

264
265
266
267
268
269
270
271
  DOFVector< double >* ProblemImplicitVec::getSignedDistance(int im , int m) 
  { 
    if ( m >= r.size() || im >= r[m].size())
	return NULL;
      return (r[m])[im]; 
    }

  DOFVector< double >* ProblemImplicitVec::getPhi1(int im, int m)
272
  {
273
274
275
276
    if (m >= phi1.size() || im >= phi1[m].size())
      return NULL;

    return (phi1[m])[im];
277
278
  }

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

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

  DOFVector< double >* ProblemImplicitVec::getLevelset(int im, int m)
  {
    if (m >= levelSet.size() || im >= levelSet[m].size())
      return NULL;

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

295
296
297
298
299
300
301
  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());
302
    for (int i = 0; i < meshes.size(); ++i)
303
304
      createImplicitMesh(i);
    return true;
305
306
  }

307

308
309
  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
310
311
312
    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
    //INFO(5,1)("initing mesh %d\n", p);
    if (!isImplicitMesh(p))
313
      return false;
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
    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);
339
340
341
    if (dofFilename.length() == 0)
      return false;

342
    double eps = getEpsilon(path, implMesh);
343
344
    if (eps < 0.0)
      return false;
345
    int serType = getType(path, implMesh);
346
347
    if (serType < 0)
      return false;
348
349
350
351

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

369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
  void ProblemImplicitVec::createMesh()
  {
    FUNCNAME("ProblemImplicitVec::createMesh()");
    ProblemVec::createMesh();
    implMesh.resize(meshes.size());
    for (int i = 0 ; i < meshes.size() ; ++i )
    {
      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;
    }

  }
396

397
398
  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, 
		                      Flag adoptFlag)
399
  {
400
     FUNCNAME("ProblemImplicitVec::initialize()");
401
     ProblemVec::initialize(initFlag);    
402
     createImplicitMesh();
403
404
  }
}