ProblemImplicit.cc 11.9 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

    GET_PARAMETER(0, path + "->macro file name", &meshFilename);
    std::string serFilename("");
    GET_PARAMETER(0, path + "->serialization file name", &serFilename);
184
    if ( meshFilename.length() > 0) {
185
      mesh->setName(path);
186
187
      readImplMesh = true;
    }
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
    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() 
  {
211
212
    FUNCNAME("ProblemImplicitScal::createImplicitMesh()");
    INFO(5,1)("setting implicit mesh\n");
213
    if (!isImplicitMesh())
214
      return false;
215
216
217
    std::string path = name + "->implicit mesh->";
    int nImplMesh(1);
    GET_PARAMETER(0, name + "nr meshes", "%d", &nImplMesh);
218

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

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

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

267

268
269
270
271
272
273
274
275
  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)
276
  {
277
278
279
280
    if (m >= phi1.size() || im >= phi1[m].size())
      return NULL;

    return (phi1[m])[im];
281
282
  }

283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
  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];
  }
298

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

311

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

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

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

373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
  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;
    }

  }
400

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