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

namespace AMDiS {
  namespace calcspace {
5
    void readDofvec(std::istream& in, DOFVector<double>* vec, Mesh* mesh) 
6
7
8
9
    {
      long size = mesh->getNumberOfVertices();
      std::string buf;
      getline(in, buf);
10
11
12
      while (buf != "vertex values: solution") 
	getline(in, buf);
      for (long i = 0; i < size; i++) {
13
14
15
16
17
        in >> buf;
        (*vec)[i] = atof(buf.c_str());
      }
    }

18

19
    void readR(std::istream& inStream, double eps, 
20
21
22
	       DOFVector<double>* r, DOFVector<double>* phi1, 
	       DOFVector<double>* phi2, DOFVector<double>* levelSet, 
	       Mesh* mesh) 
23
24
    {
      readDofvec(inStream, r, mesh);
25
26
27
28
29
30
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();
      bool checkSize = 
	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
31
      TEST_EXIT(checkSize)("something went wrong\n");
32
      for (int i = vecR.size() - 1; i >= 0; --i) {
33
34
35
36
37
38
        vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi2[i] = Phi2(vecR[i], eps);
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

39

40
    void readPhi1(istream& inStream, double eps, 
41
42
43
		  DOFVector<double>* r, DOFVector<double>* phi1, 
		  DOFVector<double>* phi2, DOFVector<double>* levelSet, 
		  Mesh* mesh) 
44
45
    {
      readDofvec(inStream, phi1, mesh);
46
47
48
49
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();
50
      
51
52
      bool checkSize = 
	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
53
      TEST_EXIT(checkSize)("something went wrong\n");
54
      for (int i = vecR.size() - 1; i>= 0; --i) {
55
56
57
58
59
60
61
        vecR[i] = Phi1ToR(vecPhi1[i], eps);
        //vecPhi2[i] = Phi2(vecR[i], eps);
        vecPhi2[i] = 1 - vecPhi1[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

62

63
    void readPhi2(istream& inStream, double eps, 
64
65
		    DOFVector<double>* r, DOFVector<double>* phi1, 
		    DOFVector<double>* phi2, DOFVector<double>* levelSet, 
66
67
68
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi2, mesh);
69
70
71
72
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();
73

74
      bool checkSize = vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
75
      TEST_EXIT(checkSize)("something went wrong\n");
76
      for (int i = vecR.size() - 1; i >= 0; --i) {
77
78
79
80
81
82
83
84
        vecR[i] = Phi2ToR(vecPhi2[i], eps);
        //vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi1[i] = 1 - vecPhi2[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }
  }

85

86
87
88
89
90
91
92
93
94
95
96
97
98
  bool ProblemImplicitScal::createImplicitMesh() 
  {
    std::string meshFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;

    std::string dofFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

99
    double eps(-1.0);
100
    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
101
    if (eps < 0)
102
103
104
105
      return false;

    int serType(-1);
    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
106
    if (serType < 0)
107
108
109
110
111
112
113
114
115
      return false;

    TEST_EXIT(mesh != NULL)("the mesh was not created\n");

    std::ifstream inStream(meshFilename.c_str());
    mesh->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin

116
    createFeSpace();
117

118
119
120
121
    r = new DOFVector< double >(getFeSpace(), "r");
    phi1 = new DOFVector<double>(getFeSpace(), "phi1");
    phi2 = new DOFVector<double>(getFeSpace(), "phi2");
    levelSet = new DOFVector< double >(getFeSpace(), "levelSet");
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
    inStream.open(dofFilename.c_str());
    switch (serType) {
    case 0:
      calcspace::readR(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 1:
      calcspace::readPhi1(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 2:
      calcspace::readPhi2(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

140

141
142
  void ProblemImplicitScal::initialize(Flag initFlag,
		  ProblemScal* adaptProblem, Flag adoptFlag)
143
  {
144
     FUNCNAME("ProblemImplicitScal::initialize()");
145
146
147
148
149
150
     ProblemScal::initialize(CREATE_MESH);
     createImplicitMesh();
     initFlag = initFlag & ~CREATE_MESH;
     ProblemScal::initialize(initFlag);
  }

151

152
153
154
155
156
157
158
  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());
159
160
    for (int i = 0; i < meshes.size(); ++i)
      createImplicitMesh(i);    
161
162
  }

163

164
165
  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
166
167
    std::string path = 
      name + "->implicit mesh[" + boost::lexical_cast< std::string >(p) + "]->";
168
169
170
171
172
173
174
175
176
177
    std::string meshFilename("");
    GET_PARAMETER(0, path + "mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;
    std::string dofFilename("");
    GET_PARAMETER(0,path + "dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

178
    double eps(-1.0);
179
    GET_PARAMETER(0, path + "eps", "%d", &eps);
180
181
    if (eps < 0.0)
      return false;
182
183
184
185
186
187
188
189
190
    int serType(-1);
    GET_PARAMETER(0, path + "type", "%d", &serType);
    if (serType < 0)
      return false;
    TEST_EXIT(meshes[p] != NULL)("the mesh was not created\n");
    std::ifstream inStream(meshFilename.c_str());
    meshes[p]->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin
191
192
193
194
195
    createFeSpace(NULL);
    r[p] = new DOFVector< double >(getFeSpace(p), "r");
    phi1[p] = new DOFVector<double>(getFeSpace(p), "phi1");
    phi2[p] = new DOFVector<double>(getFeSpace(p), "phi2");
    levelSet[p] = new DOFVector< double >(getFeSpace(p), "levelSet");
196
    inStream.open(dofFilename.c_str());
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    switch (serType) {
    case 0:
      calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p],
		       levelSet[p], meshes[p]);
      break;
    case 1:
      calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p],
			  levelSet[p], meshes[p]);
      break;
    case 2:
      calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p],
			  levelSet[p], meshes[p]);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
212
213
214
215
216
    }
    inStream.close();
    return true;
  }

217

218
219
  void ProblemImplicitVec::initialize(Flag initFlag,
		  ProblemScal* adaptProblem, Flag adoptFlag)
220
  {
221
     FUNCNAME("ProblemImplicitVec::initialize()");
222
223
224
225
226
227
     ProblemVec::initialize(CREATE_MESH);
     createImplicitMesh();
     initFlag = initFlag & ~CREATE_MESH;
     ProblemVec::initialize(initFlag);    
  }
}