BaseProblem.hh 11.1 KB
Newer Older
1
2
using namespace AMDiS;

3
4
template<typename ProblemType, bool safeguard>
BaseProblem<ProblemType, safeguard>::BaseProblem(const std::string &name_, bool createProblem) 
5
6
7
8
  : ProblemInstatBase(name_,NULL),
    prob(NULL),
    dim(1),
    dow(1),
9
10
    nTimesteps(-1),
    oldMeshChangeIdx(0),
11
    oldTimestep(0.0)
12
13
{
  // create basic problems
14
15
16
  if (createProblem)
    prob = new ProblemType(name + "->space");
  
17
  dow = Global::getGeo(WORLD);
18
}
19
20


21
22
template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::initialize(Flag initFlag,
23
24
					  ProblemStat *adoptProblem,
					  Flag adoptFlag)
25
{ FUNCNAME("BaseProblem::initialize()");
26
27

  TEST_EXIT(prob)("Problem must be created! Add flag createProblem=true to constructor!\n");
28
  
29
30
  prob->initialize(initFlag, adoptProblem, adoptFlag);
  dim = getMesh()->getDim();
31
}
32
33


34
35
template<typename ProblemType, bool safeguard>
Flag BaseProblem<ProblemType, safeguard>::initDataFromFile(AdaptInfo *adaptInfo)
36
37
38
{ 
  FUNCNAME("BaseProblem::initDataFromFile()");
  
39
  Flag initFlag;
40
  bool readDataFromFile = false;
41
42
43
  Initfile::get(name + "->read data from file", readDataFromFile, 2);
  if (!readDataFromFile)
    return initFlag;
44
  
45
46
  std::string readFormat = "arh";
  Initfile::get(name + "->read format", readFormat, 2);
Praetorius, Simon's avatar
Praetorius, Simon committed
47
  if (readFormat != "arh" && readFormat != "dat" && readFormat != "vtu" && readFormat != "multi-vtu") {
48
49
    WARNING("You can not read data from formats other than .arh, .dat or .vtu! The .arh-format is selected.\n");
  }
50
  
51
52
  // read data and mesh from arh-files/dat-files
  MSG("read data from file...\n");
53
54
  if (readFormat == "arh") 
  {
55
56
57
58
59
    std::string filename = "";
    Initfile::get(name + "->value file", filename);
    if (filename.size() == 0) return initFlag;
    if (!file_exists(filename))
      throw(std::runtime_error("The file '" + filename + "' does not exist!"));
60
61
    
    std::vector<DOFVector<double>*> solutions;
62
    for (int i = 0; i < prob->getNumComponents(); i++)
63
      solutions.push_back(prob->getSolution()->getDOFVector(i));
64
65
    ArhReader::read(filename, prob->getMesh(), solutions);
  }
66
67
  else if (readFormat == "dat") 
  {
68
69
70
71
72
73
    bool preserveMacroFileInfo = false;
    Parameters::get(prob->getMesh()->getName() + "->preserve macroFileInfo", preserveMacroFileInfo);
    if (prob->getMesh()->getMacroFileInfo() == NULL || !preserveMacroFileInfo)
      throw(std::runtime_error("Yout have to set the flag 'mesh_name->preserve macroFileInfo' to 'true' ("+boost::lexical_cast<std::string>(preserveMacroFileInfo)+"), in order to read .dat-files"));
    
    std::string filename = "";
74
    for (int i = 0; i < prob->getNumComponents(); i++) {
75
76
77
78
79
      Parameters::get(name + "->value file["+boost::lexical_cast<std::string>(i)+"]",filename);
      if (!file_exists(filename))
	throw(std::runtime_error("The file '" + filename + "'does not exist!"));
      ValueReader::readValue(filename,prob->getMesh(),prob->getSolution()->getDOFVector(i),prob->getMesh()->getMacroFileInfo());
    }
80
  }
81
82
  else if (readFormat == "vtu") 
  {
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    std::vector<std::string> filenames;
    Initfile::get(name + "->value file", filenames);
    if (filenames.size() == 0)
      return initFlag;
    int arh_idx = -1, vtu_idx = -1;
    for (size_t i = 0; i < filenames.size(); i++) {
      if (!file_exists(filenames[i]))
	throw(std::runtime_error("The file '" + filenames[i] + "' does not exist!"));
      if (filenames[i].find(".vtu") != std::string::npos)
	vtu_idx = i;
      else if (filenames[i].find(".arh") != std::string::npos)
	arh_idx = i;
      else
	throw(std::runtime_error("The file '" + filenames[i] + "' must have the extension .arh or .vtu!"));
    }
    if (arh_idx >= 0) {
      ArhReader::read(filenames[arh_idx], prob->getMesh());
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
101
102
103
104
    if (vtu_idx >= 0) {      
      std::vector<int> readComponents;
      Initfile::get(name + "->read components", readComponents);
      if (readComponents.size() == 0)
105
	for (int i = 0; i < prob->getNumComponents(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
106
107
	  readComponents.push_back(i);
	
108
109
      std::vector<DOFVector<double>*> solutions;
      std::vector<std::string> names;
Praetorius, Simon's avatar
Praetorius, Simon committed
110
111
112
      for (size_t i = 0; i < readComponents.size(); i++) {
	solutions.push_back(prob->getSolution()->getDOFVector(readComponents[i]));
	names.push_back(prob->getComponentName(readComponents[i]));
113
114
115
116
117
      }
      VtuReader::readValue(filenames[vtu_idx], prob->getMesh(), solutions, names);
    } else
      throw(std::runtime_error("You have to specify a .vtu file!"));
  }
118
119
  else if (readFormat == "multi-vtu") 
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
120
121
    size_t numFiles = 0;
    Initfile::get(name + "->number of files", numFiles);
122
    for (size_t n = 0; n < numFiles; n++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
      std::vector<std::string> filenames;
      Initfile::get(name + "->value file[" + boost::lexical_cast<std::string>(n) + "]", filenames);
      if (filenames.size() == 0)
	return initFlag;
      int arh_idx = -1, vtu_idx = -1;
      for (size_t i = 0; i < filenames.size(); i++) {
	if (!file_exists(filenames[i]))
	  throw(std::runtime_error("The file '" + filenames[i] + "' does not exist!"));
	if (filenames[i].find(".vtu") != std::string::npos)
	  vtu_idx = i;
	else if (filenames[i].find(".arh") != std::string::npos)
	  arh_idx = i;
	else
	  throw(std::runtime_error("The file '" + filenames[i] + "' must have the extension .arh or .vtu!"));
      }
      if (arh_idx >= 0) {
	ArhReader::read(filenames[arh_idx], prob->getMesh());
      }
      if (vtu_idx >= 0) {      
	std::vector<int> readComponents;
	Initfile::get(name + "->read components[" + boost::lexical_cast<std::string>(n) + "]", readComponents);
	if (readComponents.size() == 0)
145
	  for (int i = 0; i < prob->getNumComponents(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
148
149
150
151
152
153
154
155
156
157
158
	    readComponents.push_back(i);
	  
	std::vector<DOFVector<double>*> solutions;
	std::vector<std::string> names;
	for (size_t i = 0; i < readComponents.size(); i++) {
	  solutions.push_back(prob->getSolution()->getDOFVector(readComponents[i]));
	  names.push_back(prob->getComponentName(readComponents[i]));
	}
	VtuReader::readValue(filenames[vtu_idx], prob->getMesh(), solutions, names);
      } else
	throw(std::runtime_error("You have to specify a .vtu file!"));
    }
  }
159
160
161
  else {
    throw(std::runtime_error("Parameter 'read data from file' set to 'true', but no parameter 'read format' specified!"));
  }
162
  
163
164
165
166
167
168
169
170
171
172
173
  bool readPvdFromFile = false;
  Initfile::get(name + "->read pvd from file", readPvdFromFile, 2);
  if (readPvdFromFile) {
    vector<FileWriterInterface*>& fileWriters = prob->getFileWriterList();
    for (size_t i = 0; i < fileWriters.size(); i++) {
      FileWriterTemplated<double>* fileWriter = dynamic_cast<FileWriterTemplated<double>*>(fileWriters[i]);
      if (fileWriter)
	initFileWriterFromFile(adaptInfo, *fileWriter);
    }
  }
  
174
175
176
177
  initFlag.setFlag(DATA_ADOPTED);
  initFlag.setFlag(MESH_ADOPTED);

  return initFlag;
178
}
179
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
210
template<typename ProblemType, bool safeguard>
template<typename T>
void BaseProblem<ProblemType, safeguard>::initFileWriterFromFile(AdaptInfo* adaptInfo, FileWriterTemplated<T>& fileWriter, bool keep_all)
{
  using namespace pugi;
      
  if(!boost::filesystem::exists(fileWriter.getFilename() + ".pvd"))
    return;
  
  xml_document vtu;
  if(!vtu.load_file((fileWriter.getFilename() + ".pvd").c_str()))
    throw(std::runtime_error("Could not load pvd file! Error in xml structure."));

  xml_node VTKFile = vtu.child("VTKFile");
  xml_node Collection = VTKFile.child("Collection");
  
  std::vector<std::pair<double, std::string> >& paraviewAnimationFrames = fileWriter.getParaviewAnimationFrames();

  for (xml_node DataSet = Collection.child("DataSet"); DataSet; DataSet = DataSet.next_sibling("DataSet")) {
    std::string timestepStr = DataSet.attribute("timestep").value();
    std::string fileStr = DataSet.attribute("file").value();
    double time = boost::lexical_cast<double>(timestepStr);
    if (keep_all || time <= adaptInfo->getStartTime())
      paraviewAnimationFrames.push_back(std::make_pair(time, fileStr));
  }
}
  

template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::beginIteration(AdaptInfo *adaptInfo) 
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
{ FUNCNAME("BaseProblem::beginIteration()");
	
  MSG("\n");
  MSG(("[[ <"+name+"> iteration ]]\n").c_str());
    
  // assemble some blocks only once, if some conditions are fullfilled
  #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  bool fixedMatrixCell = false;
  #else
  bool fixedMatrixCell = adaptInfo->getTimestepNumber() > 0
    && std::abs(adaptInfo->getTimestep() - oldTimestep) < DBL_TOL
    && oldMeshChangeIdx == getMesh()->getChangeIndex();
  #endif
  for(size_t i = 0; i < fixedMatrixTimestep.size(); ++i) {
    prob->setAssembleMatrixOnlyOnce(
      fixedMatrixTimestep[i].first, 
      fixedMatrixTimestep[i].second,
      fixedMatrixCell);
  }
230
}
231
232


233
234
template<typename ProblemType, bool safeguard>
Flag BaseProblem<ProblemType, safeguard>::oneIteration(AdaptInfo *adaptInfo, 
235
                                Flag toDo) 
236
237
{ 
  FUNCNAME("BaseProblem::oneIteration()");
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

  Flag flag, markFlag;
  
  if (toDo.isSet(MARK))
    markFlag = prob->markElements(adaptInfo);
  else
    markFlag = 3;

  // refine
  if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
    flag = prob->refineMesh(adaptInfo);  
  // coarsen
  if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
    flag |= prob->coarsenMesh(adaptInfo);

  #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  MPI::COMM_WORLD.Barrier();
  #endif

257
258
#ifndef NONLIN_PROBLEM
  if (toDo.isSet(BUILD))
259
    prob->buildAfterCoarsen(adaptInfo, markFlag, true, true);
260
#endif
261
262

  if (toDo.isSet(SOLVE)) {
263
    if (safeguard) {
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
      try {
	prob->solve(adaptInfo);
      } catch(std::exception &e) {
	WARNING(("Could not solve system. Simulation will be stoped here! ERROR: "
	  + std::string(e.what()) + "\n").c_str());
	adaptInfo->setTime(adaptInfo->getEndTime());
	adaptInfo->setTimestepNumber(adaptInfo->getNumberOfTimesteps());
	return flag;
      } catch(...) {
	WARNING("Could not solve system. Simulation will be stoped here! Unknown ERROR\n");
	adaptInfo->setTime(adaptInfo->getEndTime());
	adaptInfo->setTimestepNumber(adaptInfo->getNumberOfTimesteps());
	return flag;
      }
    } else {
      prob->solve(adaptInfo);
    }
  }

  if (toDo.isSet(ESTIMATE))
    prob->estimate(adaptInfo);

  oldTimestep = adaptInfo->getTimestep();
  oldMeshChangeIdx = getMesh()->getChangeIndex();

  return flag;
290
}
291
292


293
294
template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::endIteration(AdaptInfo *adaptInfo) 
295
296
297
298
{ FUNCNAME("BaseProblem::endIteration()");
	
  MSG("\n");
  MSG(("[[ end of <"+name+"> iteration ]]\n").c_str());
299
}
300
301


302
303
template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::closeTimestep(AdaptInfo *adaptInfo)
304
{
305
  writeFiles(adaptInfo, false);
306
}
307
308


309
310
template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::addTimeOperator(ProblemStat *prob, int i, int j)
311
{
312
313
314
315
316
317
318
  Operator *op = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
    op->addZeroOrderTerm(new Simple_ZOT);
  Operator *opRhs = new Operator(prob->getFeSpace(i));
    opRhs->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(j), NULL));

  prob->addMatrixOperator(op, i, j, getInvTau(), getInvTau());
  prob->addVectorOperator(opRhs, i, getInvTau(), getInvTau());
319
}
320
321


322
323
template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::addTimeOperator(RosenbrockStationary *prob, int i, int j)
324
{
325
  prob->addTimeOperator(i,j);
326
}