CouplingProblemStat.h 12.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21
22
23
24
25
26
27
28
29
30
31
32
33



/** \file CouplingProblemStat.h */

#ifndef AMDIS_COUPLING_PROBLEM_STAT_H
#define AMDIS_COUPLING_PROBLEM_STAT_H

#include <vector>
#include <set>
#include <list>
#include "AMDiS_fwd.h"
#include "ProblemStat.h"
#include "Initfile.h"
34
#include "utility/to_string.hpp"
35
36
37

namespace AMDiS {

38
  namespace detail
39
40
  {

41
42
43
44
45
    /** \brief
    * This class defines a coupled stationary problem definition in sequential
    * computations. 
    */
    template<typename ProblemStatType>
46
    class CouplingProblemStat : public ProblemStatSeq
47
    {
48
49
50
51
52
53
54
55
56
57
58
    protected:
      typedef ProblemStatSeq  super;
      
      using super::nComponents;
      using super::meshes;
      using super::nMeshes;
      using super::feSpaces;
      using super::name;
      using super::refinementManager;
      using super::coarseningManager;
      
59
60
61
    public:
      /// Constructor
      CouplingProblemStat(std::string name_)
62
63
	: super(name_),
	  dim(-1)
64
65
66
67
68
69
70
71
72
73
74
      {}

      /// add problem by number
      virtual void addProblem(ProblemStatType* prob)
      {
	problems.push_back(prob);
	nComponents += prob->getNumComponents();
      };

      /// Initialisation of the problem.
      virtual void initialize(Flag initFlag,
75
			      ProblemStatSeq *adoptProblem = NULL,
76
			      Flag adoptFlag = INIT_NOTHING)
77
      { FUNCNAME("CouplingProblemStat::initialize()");
78

79
80
81
82
83
	super::initialize(initFlag - INIT_MESH);
	
	const Flag DEFAULT_INIT = (INIT_FE_SPACE | INIT_MESH | CREATE_MESH | INIT_SYSTEM | INIT_SOLVER | INIT_ESTIMATOR | INIT_MARKER | INIT_FILEWRITER);
	for (size_t p = 0; p < problems.size(); ++p) {
	  problems[p]->initialize(initFlag - DEFAULT_INIT);
84
85
	}

86
87
88
89
90
91
92
93
94
95
	for (size_t i = 0; i < meshes.size(); i++) {
	  int globalRefinements = 0;

	  // If AMDiS is compiled for parallel computations, the global refinements are
	  // ignored here. Later, each rank will add the global refinements to its
	  // private mesh.
  #ifndef HAVE_PARALLEL_DOMAIN_AMDIS
	  Parameters::get(meshes[i]->getName() + "->global refinements",
			  globalRefinements);
  #endif
96

97
	  bool initMesh = initFlag.isSet(INIT_MESH);
98

99
100
101
102
103
104
105
106
107
108
109
110
	  // Initialize the meshes if there is no serialization file.
	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized())) {
	    meshes[i]->initialize();
	    refinementManager->globalRefine(meshes[i], globalRefinements);
	  }
	}
      }
      
      
      /// Used in \ref initialize().
      virtual void createMesh() override
      {
111
	// all problems must have the same dimension (?)
112
	dim = 0;
113
114
115
116
	Parameters::get(name + "->dim", dim);
	TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
		      name.c_str());

117
118
119
	std::map<std::pair<std::string, int>, Mesh*> meshByName;  // (name, refSet) --> Mesh*
	typedef std::map<std::pair<std::string, int>, Mesh*>::iterator MeshIterator;

120
	std::vector<std::set<Mesh*> > meshesForProblems(problems.size());
121
	
122
123
	for (size_t i = 0; i < problems.size(); ++i) {
	  TEST_EXIT(problems[i])("problem[%d] does not exist!\n",i);
124
125
126
127
128
129
130
131
132
133
	  
	  int nComponents = problems[i]->getNumComponents();
	  
	  int nAddComponents = 0;
	  Parameters::get(problems[i]->getName() + "->additional components", nAddComponents);
	  
	  problems[i]->componentMeshes.resize(nComponents + nAddComponents);
	
	  for (size_t j = 0; j < nComponents + nAddComponents; j++) {
	    // name of the mesh
134
135
136
137
	    std::string meshName("");
	    Parameters::get(problems[i]->getName() + "->mesh", meshName);
	    TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n",
				      problems[i]->getName().c_str());
138
139
140
141
142
143
144
145
146
147
148
149
150
151
				      
	    // dimension of the mesh
	    int mesh_dim = 0;
	    Parameters::get(problems[i]->getName() + "->dim", mesh_dim);
	    TEST_EXIT(dim == mesh_dim)("Mesh-dimension must be the same for all problems!\n");
		
	    // refinement set (optional)
	    int refSet = 0;
	    Parameters::get(problems[i]->getName() + "->refinement set[" + to_string(j) + "]", refSet);

	    // create a new Mesh only if not already created for other problem
	    Mesh* componentMesh;
	    MeshIterator meshIt = meshByName.find(std::make_pair(meshName, refSet));
	    if (meshIt == meshByName.end()) {
152
153
	      Mesh *newMesh = new Mesh(meshName, dim);
	      meshes.push_back(newMesh);
154
155
	      meshByName[std::make_pair(meshName, refSet)] = newMesh;	    
	      componentMesh = newMesh;
156
	      nMeshes++;
157
158
159
160
161
162
163
164
165
166
167
168
169
170
	    } else {
	      componentMesh = meshIt->second;
	    }
	    problems[i]->componentMeshes[j] = componentMesh;
	  }
	  
	  // copy unqiue set of meshes to problem[i]->meshes
	  std::set<Mesh*> uniqueMeshes;
	  for (size_t j = 0; j < problems[i]->componentMeshes.size(); ++j)
	    uniqueMeshes.insert(problems[i]->componentMeshes[j]);
	  problems[i]->meshes.clear();
	  problems[i]->meshes.insert(problems[i]->meshes.begin(), uniqueMeshes.begin(), uniqueMeshes.end());
	}
      }
171

172
173
174
175
176
177
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
211
212
213
214
215
216
217
218
219
220
221
      /// Used in \ref initialize().
      virtual void createFeSpace(DOFAdmin *admin) override
      {
	std::vector<std::set<FiniteElemSpace const*> > feSpacesForProblems(problems.size());
	std::map<std::pair<Mesh*, std::string>, FiniteElemSpace*> feSpaceMap;
	
	for (size_t p = 0; p < problems.size(); ++p) {
	  TEST_EXIT(problems[p])("problem[%d] does not exist!\n",p);
	  
	  int nComponents = problems[p]->getNumComponents();
	  
	  int nAddComponents = 0;
	  Parameters::get(problems[p]->getName() + "->additional components", nAddComponents);
	  problems[p]->componentSpaces.resize(nComponents + nAddComponents, NULL);
	  problems[p]->traverseInfo.resize(nComponents);
	
	  for (size_t i = 0; i < nComponents + nAddComponents; i++) {
	  
	    std::string componentString = "[" + to_string(i) + "]";
	    
	    std::string feSpaceName = "";
	    std::string initFileStr = problems[p]->getName() + "->feSpace" + componentString;
	    Parameters::get(initFileStr, feSpaceName);
	    
	    // synonym for "feSpace"
	    if (feSpaceName.size() == 0) {
	      initFileStr = problems[p]->getName() + "->finite element space" + componentString;
	      Parameters::get(initFileStr, feSpaceName);
	    }
	    
	    // for backward compatibility also accept the old syntax
	    if (feSpaceName.size() == 0) {
	      int degree = 1;
	      initFileStr = problems[p]->getName() + "->polynomial degree" + componentString;
	      Parameters::get(initFileStr, degree);
	      TEST_EXIT(degree > 0)
		("Poynomial degree in component %d must be larger than zero!\n", i);
		
	      feSpaceName = "Lagrange" + to_string(degree);
	    } 
	    
	    if (feSpaceName.size() == 0)
	      feSpaceName = "Lagrange1";

	    if (feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)] == NULL) {      
	      BasisFunctionCreator *basisFctCreator = 
		dynamic_cast<BasisFunctionCreator*>(CreatorMap<BasisFunction>::getCreator(feSpaceName, initFileStr));
	      TEST_EXIT(basisFctCreator)
		("No valid basisfunction type found in parameter \"%s\"\n", initFileStr.c_str());
	      basisFctCreator->setDim(dim);
222
223
	      
	      FiniteElemSpace *newFeSpace = 
224
225
226
227
228
		FiniteElemSpace::provideFeSpace(admin, basisFctCreator->create(),
						problems[p]->componentMeshes[i], 
						"FeSpace" + componentString + " (" + feSpaceName + ")");

	      feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)] = newFeSpace;
229
230
231
	      feSpaces.push_back(newFeSpace);
	    }

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
	    problems[p]->componentSpaces[i] = feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)];
	  }
	  
	  // copy unqiue set of meshes to problem[p]->meshes
	  std::set<FiniteElemSpace const*> uniqueFeSpaces;
	  for (size_t i = 0; i < problems[p]->componentSpaces.size(); ++i)
	    uniqueFeSpaces.insert(problems[p]->componentSpaces[i]);
	  problems[p]->feSpaces.clear();
	  problems[p]->feSpaces.insert(problems[p]->feSpaces.begin(), uniqueFeSpaces.begin(), uniqueFeSpaces.end());
	  
	  // create traverseInfo
	  for (int i = 0; i < nComponents; i++) {
	    for (int j = 0; j < nComponents; j++)
	      problems[p]->traverseInfo.getMatrix(i, j).setFeSpace(problems[p]->componentSpaces[i], problems[p]->componentSpaces[j]);
	    
	    problems[p]->traverseInfo.getVector(i).setFeSpace(problems[p]->componentSpaces[i]);
248
249
	  }
	}
250

251
	// create dof admin for vertex dofs if neccessary
252
	for (size_t i = 0; i < meshes.size(); i++) {
253
254
255
256
	  if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
	    DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
	    ln_dof[VERTEX] = 1;
	    meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
257
	  }
258
	}
259
260
      }

261
262
      /// Used in \ref initialize().
      virtual void createRefCoarseManager() override
263
      {
264
	FUNCNAME("CouplingProblemStat::createRefCoarseManager()");
265
	assert( dim > 0);
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282

	switch (dim) {
	case 1:
	  coarseningManager = new CoarseningManager1d();
	  refinementManager = new RefinementManager1d();
	  break;
	case 2:
	  coarseningManager = new CoarseningManager2d();
	  refinementManager = new RefinementManager2d();
	  break;
	case 3:
	  coarseningManager = new CoarseningManager3d();
	  refinementManager = new RefinementManager3d();
	  break;
	default:
	  ERROR_EXIT("invalid dim!\n");
	}
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
	
	for (size_t p = 0; p < problems.size(); p++) {
	  problems[p]->setRefinementManager(refinementManager);
	  problems[p]->setCoarseningManager(coarseningManager);
	}
      }

      /// Used in \ref initialize().
      virtual void createMatricesAndVectors() override
      {
	for (size_t p = 0; p < problems.size(); ++p) {
	  assert( problems[p] );
	  problems[p]->createMatricesAndVectors();
	}
      }

      /// Used in \ref initialize().
      virtual void createSolver() override
      {
	for (size_t p = 0; p < problems.size(); ++p) {
	  assert( problems[p] );
	  problems[p]->createSolver();
	}
306
307
      }

308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
      /// Used in \ref initialize().
      virtual void createEstimator() override
      {
	for (size_t p = 0; p < problems.size(); ++p) {
	  assert( problems[p] );
	  problems[p]->createEstimator();
	}
      }

      /// Used in \ref initialize().
      virtual void createMarker() override
      {
	for (size_t p = 0; p < problems.size(); ++p) {
	  assert( problems[p] );
	  problems[p]->createMarker();
	}
      }

      /// Used in \ref initialize().
      virtual void createFileWriter() override
      {
	for (size_t p = 0; p < problems.size(); ++p) {
	  assert( problems[p] );
	  problems[p]->createFileWriter();
	}
      }


336
337
338
339
340
      /// Returns number of managed problems
      virtual int getNumProblems() 
      { 
	return problems.size(); 
      }
341

342
343
344
345
346
347
348
349
      /** \brief
      * Returns the problem with the given number. If only one problem
      * is managed by this master problem, the number hasn't to be given.
      */
      virtual ProblemStatType *getProblem(int number = 0)
      { 
	return problems[number]; 
      }
350

351
352
353
354
355
      /// Returns \ref meshes[i]
      inline Mesh* getMesh(int number = 0) 
      {
	return meshes[number]; 
      }
356

357
358
359
360
      using super::getNumComponents;
      using super::getRefinementManager;
      using super::getCoarseningManager;
      using super::getName;
361

362
363
    protected:
      
364
365
366
367
368
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
396
397
398
//       /// Name of this problem.
//       std::string name;
// 
//       /// Number of problem components
//       int nComponents;
      
      /// unqiue mesh-dimension for all problems
      int dim;

//       /** \brief
//       * Number of problem meshes. If all components are defined on the same mesh,
//       * this number is 1. Otherwise, this variable is the number of different meshes
//       * within the problem.
//       */
//       int nMeshes;

//       /// FE spaces of this problem.
//       std::vector<FiniteElemSpace*> feSpaces;
// 
//       /// Meshes of this problem.
//       std::vector<Mesh*> meshes;

//       /** \brief
//       * All actions of mesh refinement are performed by refinementManager.
//       * If new refinement algorithms should be realized, one has to override
//       * RefinementManager and give one instance of it to AdaptStationary.
//       */
//       RefinementManager *refinementManager;
// 
//       /** \brief
//       * All actions of mesh coarsening are performed by coarseningManager.
//       * If new coarsening algorithms should be realized, one has to override
//       * CoarseningManager and give one instance of it to AdaptStationary.
//       */
//       CoarseningManager *coarseningManager;
399
400
401

      std::vector<ProblemStatType*> problems;
    };
402

403
  } // end namespace detail
404

405
406
407
  typedef detail::CouplingProblemStat<ProblemStat> CouplingProblemStat;
  
} // end namespace AMDiS
408
409

#endif