CouplingProblemStat.h 11.4 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * 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.
18
 *
19
 ******************************************************************************/
20 21 22 23 24 25 26 27 28 29 30 31 32



/** \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"
Praetorius, Simon's avatar
Praetorius, Simon committed
33 34
#include "CoarseningManager.h"
#include "RefinementManager.h"
35
#include "Initfile.h"
36
#include "utility/to_string.hpp"
37 38 39

namespace AMDiS {

40
  namespace detail
41 42
  {

43 44
    /** \brief
    * This class defines a coupled stationary problem definition in sequential
45
    * computations.
46 47
    */
    template<typename ProblemStatType>
48
    class CouplingProblemStat : public ProblemStatSeq
49
    {
50 51
    protected:
      typedef ProblemStatSeq  super;
52

53 54 55 56 57 58 59
      using super::nComponents;
      using super::meshes;
      using super::nMeshes;
      using super::feSpaces;
      using super::name;
      using super::refinementManager;
      using super::coarseningManager;
60

61 62 63
    public:
      /// Constructor
      CouplingProblemStat(std::string name_)
64 65
	: super(name_),
	  dim(-1)
66 67 68 69 70 71 72 73 74 75
      {}

      /// add problem by number
      virtual void addProblem(ProblemStatType* prob)
      {
	problems.push_back(prob);
      };

      /// Initialisation of the problem.
      virtual void initialize(Flag initFlag,
76
			      ProblemStatSeq *adoptProblem = NULL,
77
			      Flag adoptFlag = INIT_NOTHING)
78
      {
79
	super::initialize(initFlag - INIT_MESH);
80

81 82 83
	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
	  // Initialize the meshes if there is no serialization file.
	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized())) {
	    meshes[i]->initialize();
	    refinementManager->globalRefine(meshes[i], globalRefinements);
	  }
	}
      }
106 107


108 109 110
      /// 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
	  int nComponents = problems[i]->getNumComponents();
126

127 128
	  int nAddComponents = 0;
	  Parameters::get(problems[i]->getName() + "->additional components", nAddComponents);
129

130
	  problems[i]->componentMeshes.resize(nComponents + nAddComponents);
131

132
	  for (int j = 0; j < nComponents + nAddComponents; j++) {
133
	    // 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
	    // 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");
143

144 145 146 147 148 149 150 151
	    // 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
	      meshByName[std::make_pair(meshName, refSet)] = newMesh;
155
	      componentMesh = newMesh;
156
	      nMeshes++;
157 158 159 160 161
	    } else {
	      componentMesh = meshIt->second;
	    }
	    problems[i]->componentMeshes[j] = componentMesh;
	  }
162

163 164 165 166 167 168 169 170
	  // 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
      /// 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;
177

178 179
	for (size_t p = 0; p < problems.size(); ++p) {
	  TEST_EXIT(problems[p])("problem[%d] does not exist!\n",p);
180

181
	  int nComponents = problems[p]->getNumComponents();
182

183 184 185 186
	  int nAddComponents = 0;
	  Parameters::get(problems[p]->getName() + "->additional components", nAddComponents);
	  problems[p]->componentSpaces.resize(nComponents + nAddComponents, NULL);
	  problems[p]->traverseInfo.resize(nComponents);
187

188
	  for (int i = 0; i < nComponents + nAddComponents; i++) {
189

190
	    std::string componentString = "[" + to_string(i) + "]";
191

192 193 194
	    std::string feSpaceName = "";
	    std::string initFileStr = problems[p]->getName() + "->feSpace" + componentString;
	    Parameters::get(initFileStr, feSpaceName);
195

196 197 198 199 200
	    // synonym for "feSpace"
	    if (feSpaceName.size() == 0) {
	      initFileStr = problems[p]->getName() + "->finite element space" + componentString;
	      Parameters::get(initFileStr, feSpaceName);
	    }
201

202 203 204 205 206 207 208
	    // 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);
209

210
	      feSpaceName = "Lagrange" + to_string(degree);
211 212
	    }

213 214 215
	    if (feSpaceName.size() == 0)
	      feSpaceName = "Lagrange1";

216 217
	    if (feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)] == NULL) {
	      BasisFunctionCreator *basisFctCreator =
218 219 220 221
		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
		FiniteElemSpace::provideFeSpace(admin, basisFctCreator->create(),
225
						problems[p]->componentMeshes[i],
226 227 228
						"FeSpace" + componentString + " (" + feSpaceName + ")");

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

232 233
	    problems[p]->componentSpaces[i] = feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)];
	  }
234

235 236 237 238 239 240
	  // 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());
241

242 243 244 245
	  // 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]);
246

247
	    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
	  if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
	    DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
	    ln_dof[VERTEX] = 1;
256
	    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
      /// Returns number of managed problems
337 338 339
      virtual int getNumProblems()
      {
	return problems.size();
340
      }
341

342 343 344 345 346
      /** \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)
347 348
      {
	return problems[number];
349
      }
350

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

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

362
    protected:
363 364 365
      /// unqiue mesh-dimension for all problems
      int dim;

366 367
      std::vector<ProblemStatType*> problems;
    };
368

369
  } // end namespace detail
370

371
  typedef detail::CouplingProblemStat<ProblemStat> CouplingProblemStat;
372

373
} // end namespace AMDiS
374 375

#endif