ProblemImplicit.cc 7.12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include "ProblemImplicit.h"

namespace AMDiS {
  namespace calcspace {
    void readDofvec(std::istream& in, DOFVector< double >* vec, Mesh* mesh) 
    {
      long size = mesh->getNumberOfVertices();
      std::string buf;
      getline(in, buf);
      while(buf != "vertex values: solution") getline(in, buf);
      for(long i = 0; i < size; i++)
      {
        in >> buf;
        (*vec)[i] = atof(buf.c_str());
      }
    }

    void readR(std::istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, r, mesh);
      std::vector< double >& vecR = r->getVector();
      std::vector< double >& vecPhi1 = phi1->getVector();
      std::vector< double >& vecPhi2 = phi2->getVector();
      std::vector< double >& vecLevelSet = levelSet->getVector();
28
29
30
      bool checkSize = vecR.size() != vecPhi1.size() &&
	      vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
      for (int i = vecR.size()-1; i >= 0; --i) 
      {
        vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi2[i] = Phi2(vecR[i], eps);
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

    void readPhi1(istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi1, mesh);
      std::vector< double >& vecR = r->getVector();
      std::vector< double >& vecPhi1 = phi1->getVector();
      std::vector< double >& vecPhi2 = phi2->getVector();
      std::vector< double >& vecLevelSet = levelSet->getVector();
49
50
51
52
      
      bool checkSize = vecR.size() != vecPhi1.size() &&
	      vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
      for (int i = vecR.size()-1; i>= 0; --i) 
      {
        vecR[i] = Phi1ToR(vecPhi1[i], eps);
        //vecPhi2[i] = Phi2(vecR[i], eps);
        vecPhi2[i] = 1 - vecPhi1[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

    void readPhi2(istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi2, mesh);
      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
75
      bool checkSize = vecR.size() != vecPhi1.size() &&
	      vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
76
77
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
105
106
107
108
109
110
111
112
113
114
115
116
117
      for (int i = vecR.size()-1; i >= 0; --i) 
      {
        vecR[i] = Phi2ToR(vecPhi2[i], eps);
        //vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi1[i] = 1 - vecPhi2[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }


  }

  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;

    double eps(-1);
    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
    if (eps<0)
      return false;

    int serType(-1);
    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
    if (serType<0)
      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

118
    createFeSpace();
119

120
121
122
123
    r = new DOFVector< double >(getFeSpace(), "r");
    phi1 = new DOFVector<double>(getFeSpace(), "phi1");
    phi2 = new DOFVector<double>(getFeSpace(), "phi2");
    levelSet = new DOFVector< double >(getFeSpace(), "levelSet");
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    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;
  }

142
143
  void ProblemImplicitScal::initialize(Flag initFlag,
		  ProblemScal* adaptProblem, Flag adoptFlag)
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
  {
     FUNCNAME("ImplicitScal::initialize()");
     ProblemScal::initialize(CREATE_MESH);
     createImplicitMesh();
     initFlag = initFlag & ~CREATE_MESH;
     ProblemScal::initialize(initFlag);
  }

  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());
    for(int i=0;i<meshes.size();++i) {
      createImplicitMesh(i);
    }
  }

  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
178
179
180
181
182
183
184
185
186
187
188
189
190
    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;

    double eps(-1);
    GET_PARAMETER(0, path + "eps", "%d", &eps);
    if(eps < 0)
        return false;
    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
197
198
199
    inStream.open(dofFilename.c_str());
    switch (serType) 
    {
      case 0:
200
201
        calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p],
		       	levelSet[p], meshes[p]);
202
203
        break;
      case 1:
204
205
	calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p],
		       	levelSet[p], meshes[p]);
206
207
        break;
      case 2:
208
209
	calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p],
		       	levelSet[p], meshes[p]);
210
211
212
213
214
215
216
217
        break;
      default:
        WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

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