diff --git a/extensions/demo/src/drivenCavity.cc b/extensions/demo/src/drivenCavity.cc
index e1a622a379d56f4952677673ed5442bef0f1351c..00b76a5362a336eedfe724f5b4901d877912a97e 100644
--- a/extensions/demo/src/drivenCavity.cc
+++ b/extensions/demo/src/drivenCavity.cc
@@ -13,8 +13,16 @@ using namespace boost::posix_time;
struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
{
double operator()(const WorldVector<double> &x) const {
- return std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
+ double vel = std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
+
+ return (*time < 0.1 ? sin((*time) * m_pi/0.2)*vel : vel);
}
+ void setTimePtr(double* time_)
+ {
+ time = time_;
+ }
+private:
+ double* time;
};
class CHNS_DrivenCavity : public CahnHilliardNavierStokes
@@ -40,6 +48,8 @@ public:
transformDOF(prob->getSolution()->getDOFVector(0),
new SignedDistToCh(initialEps));
}
+
+ drivenCavityBC->setTimePtr(adaptInfo->getTimePtr());
}
void fillBoundaryConditions()
@@ -53,10 +63,14 @@ public:
for (size_t i = 0; i < dow; i++)
getProblem(0)->addDirichletBC(1, 2+i, 2+i, zeroDOF);
+ drivenCavityBC = new DrivenCavityBC;
/// at upper wall: prescribed velocity
- getProblem(0)->addDirichletBC(2, 2, 2, new DrivenCavityBC);
+ getProblem(0)->addDirichletBC(2, 2, 2, drivenCavityBC);
getProblem(0)->addDirichletBC(2, 3, 3, zeroDOF);
}
+
+private:
+ DrivenCavityBC* drivenCavityBC;
};
@@ -96,7 +110,7 @@ public:
/// Set initial condition and perform initial refinement
virtual void solveInitialProblem(AdaptInfo *adaptInfo)
{
- refFunction= new PhaseFieldRefinement("mesh->refinement",chnsProb->getMesh());
+ refFunction= new PhaseFieldRefinement(chnsProb->getMesh());
refinement= new RefinementLevelDOF(
chnsProb->getProblem(0)->getFeSpace(0),
refFunction,
diff --git a/extensions/demo/src/drivenCavity_rb.cc b/extensions/demo/src/drivenCavity_rb.cc
new file mode 100644
index 0000000000000000000000000000000000000000..91b41dbaf7555310177fdd33cdc358684d33fdbe
--- /dev/null
+++ b/extensions/demo/src/drivenCavity_rb.cc
@@ -0,0 +1,161 @@
+#include "AMDiS.h"
+#include "CahnHilliardNavierStokes_RB.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#include "time/ExtendedRosenbrockAdaptInstationary.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
+{
+ double operator()(const WorldVector<double> &x) const {
+ double vel = std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
+
+ return vel;
+ }
+};
+
+class CHNS_DrivenCavity : public CahnHilliardNavierStokes_RB
+{
+public:
+ CHNS_DrivenCavity(std::string name_) : CahnHilliardNavierStokes_RB(name_) {}
+
+ void solveInitialProblem(AdaptInfo* adaptInfo)
+ {
+ super::solveInitialProblem(adaptInfo);
+ /// horizontale Linie
+ double a= 0.0, dir= -1.0;
+ double initialEps = eps;
+ Initfile::get(name + "->initial epsilon", initialEps);
+ Initfile::get(name + "->line->pos", a);
+ Initfile::get(name + "->line->direction", dir);
+
+ /// create phase-field from signed-dist-function
+ if (doubleWell == 0) {
+ prob->getSolution()->getDOFVector(0)->interpol(new SignedDistFctToPhaseField(initialEps, new Plane(a, dir)));
+ } else {
+ prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
+ transformDOF(prob->getSolution()->getDOFVector(0),
+ new SignedDistToCh(initialEps));
+ }
+ }
+
+ void fillBoundaryConditions()
+ { FUNCNAME("NS_DrivenCavity::fillBoundaryConditions()");
+
+ AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+ size_t dow = Global::getGeo(WORLD);
+
+ /// at rigid wall: no-slip boundary condition
+ for (size_t i = 0; i < dow; i++)
+ getProblem(0)->addDirichletBC(1, 2+i, 2+i, zero);
+
+ drivenCavityBC = new DrivenCavityBC;
+ /// at upper wall: prescribed velocity
+ getProblem(0)->addDirichletBC(2, 2, 2, drivenCavityBC);
+ getProblem(0)->addDirichletBC(2, 3, 3, zero);
+ }
+
+private:
+ DrivenCavityBC* drivenCavityBC;
+};
+
+
+class RefinementTimeInterface : public CouplingTimeInterface
+{
+public:
+
+ RefinementTimeInterface(CHNS_DrivenCavity *chnsProb_) :
+ chnsProb(chnsProb_),
+ refFunction(NULL),
+ refinement(NULL)
+ {
+ addTimeInterface(chnsProb);
+ }
+
+ ~RefinementTimeInterface()
+ {
+ if (refFunction)
+ delete refFunction;
+ if (refinement)
+ delete refinement;
+ }
+
+ virtual void initTimeInterface()
+ {
+ chnsProb->initTimeInterface();
+ }
+
+ /// Called at the end of each timestep.
+ virtual void closeTimestep(AdaptInfo *adaptInfo)
+ {
+ CouplingTimeInterface::closeTimestep(adaptInfo);
+
+ refinement->refine(1);
+ }
+
+ /// Set initial condition and perform initial refinement
+ virtual void solveInitialProblem(AdaptInfo *adaptInfo)
+ {
+ refFunction= new PhaseFieldRefinement(chnsProb->getMesh());
+ refinement= new RefinementLevelDOF(
+ chnsProb->getProblem(0)->getFeSpace(0),
+ refFunction,
+ new PhaseDOFView<double>(chnsProb->getProblem(0)->getSolution()->getDOFVector(0)));
+
+ // initial refinement
+ refinement->refine(0);
+
+ bool initialRefinement = true;
+ Parameters::get(chnsProb->getName() + "->initial refinement", initialRefinement);
+
+ if (initialRefinement) {
+ // refine until interfaces is solved
+ for (int i = 0; i < 3; ++i) {
+ chnsProb->solveInitialProblem(adaptInfo);
+ refinement->refine((i < 4 ? 4 : 10));
+ }
+ }
+
+ CouplingTimeInterface::solveInitialProblem(adaptInfo);
+ }
+
+protected:
+
+ CHNS_DrivenCavity *chnsProb;
+ PhaseFieldRefinement* refFunction;
+ RefinementLevelDOF *refinement;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+ AMDiS::init(argc, argv);
+
+ CHNS_DrivenCavity chnsProb("chns");
+ chnsProb.initialize(INIT_ALL);
+
+ RefinementTimeInterface timeInterface(&chnsProb);
+
+ // Adapt-Infos
+ AdaptInfo adaptInfo("adapt", chnsProb.getNumComponents());
+ ExtendedRosenbrockAdaptInstationary<CHNS_DrivenCavity> adaptInstat("adapt", chnsProb, adaptInfo, timeInterface, adaptInfo);
+
+
+ ptime start_time = microsec_clock::local_time();
+ chnsProb.initTimeInterface();
+ int error_code = adaptInstat.adapt();
+ time_duration td = microsec_clock::local_time()-start_time;
+
+ MSG("elapsed time= %d sec\n", td.total_seconds());
+
+ AMDiS::finalize();
+
+ return error_code;
+};
diff --git a/extensions/demo/src/drivenCavity_twophase_rb.cc b/extensions/demo/src/drivenCavity_twophase_rb.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3c1ee313954ab6655011dee10318483d1e7cbf30
--- /dev/null
+++ b/extensions/demo/src/drivenCavity_twophase_rb.cc
@@ -0,0 +1,161 @@
+#include "AMDiS.h"
+#include "CahnHilliardNavierStokes_TwoPhase_RB.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#include "time/ExtendedRosenbrockAdaptInstationary.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
+{
+ double operator()(const WorldVector<double> &x) const {
+ double vel = std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
+
+ return vel;
+ }
+};
+
+class CHNS_DrivenCavity : public CahnHilliardNavierStokes_TwoPhase_RB
+{
+public:
+ CHNS_DrivenCavity(std::string name_) : CahnHilliardNavierStokes_TwoPhase_RB(name_) {}
+
+ void solveInitialProblem(AdaptInfo* adaptInfo)
+ {
+ super::solveInitialProblem(adaptInfo);
+ /// horizontale Linie
+ double a= 0.0, dir= -1.0;
+ double initialEps = eps;
+ Initfile::get(name + "->initial epsilon", initialEps);
+ Initfile::get(name + "->line->pos", a);
+ Initfile::get(name + "->line->direction", dir);
+
+ /// create phase-field from signed-dist-function
+ if (doubleWell == 0) {
+ prob->getSolution()->getDOFVector(0)->interpol(new SignedDistFctToPhaseField(initialEps, new Plane(a, dir)));
+ } else {
+ prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
+ transformDOF(prob->getSolution()->getDOFVector(0),
+ new SignedDistToCh(initialEps));
+ }
+ }
+
+ void fillBoundaryConditions()
+ { FUNCNAME("NS_DrivenCavity::fillBoundaryConditions()");
+
+ AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+ size_t dow = Global::getGeo(WORLD);
+
+ /// at rigid wall: no-slip boundary condition
+ for (size_t i = 0; i < dow; i++)
+ getProblem(0)->addDirichletBC(1, 2+i, 2+i, zero);
+
+ drivenCavityBC = new DrivenCavityBC;
+ /// at upper wall: prescribed velocity
+ getProblem(0)->addDirichletBC(2, 2, 2, drivenCavityBC);
+ getProblem(0)->addDirichletBC(2, 3, 3, zero);
+ }
+
+private:
+ DrivenCavityBC* drivenCavityBC;
+};
+
+
+class RefinementTimeInterface : public CouplingTimeInterface
+{
+public:
+
+ RefinementTimeInterface(CHNS_DrivenCavity *chnsProb_) :
+ chnsProb(chnsProb_),
+ refFunction(NULL),
+ refinement(NULL)
+ {
+ addTimeInterface(chnsProb);
+ }
+
+ ~RefinementTimeInterface()
+ {
+ if (refFunction)
+ delete refFunction;
+ if (refinement)
+ delete refinement;
+ }
+
+ virtual void initTimeInterface()
+ {
+ chnsProb->initTimeInterface();
+ }
+
+ /// Called at the end of each timestep.
+ virtual void closeTimestep(AdaptInfo *adaptInfo)
+ {
+ CouplingTimeInterface::closeTimestep(adaptInfo);
+
+ refinement->refine(1);
+ }
+
+ /// Set initial condition and perform initial refinement
+ virtual void solveInitialProblem(AdaptInfo *adaptInfo)
+ {
+ refFunction= new PhaseFieldRefinement(chnsProb->getMesh());
+ refinement= new RefinementLevelDOF(
+ chnsProb->getProblem(0)->getFeSpace(0),
+ refFunction,
+ new PhaseDOFView<double>(chnsProb->getProblem(0)->getSolution()->getDOFVector(0)));
+
+ // initial refinement
+ refinement->refine(0);
+
+ bool initialRefinement = true;
+ Parameters::get(chnsProb->getName() + "->initial refinement", initialRefinement);
+
+ if (initialRefinement) {
+ // refine until interfaces is solved
+ for (int i = 0; i < 3; ++i) {
+ chnsProb->solveInitialProblem(adaptInfo);
+ refinement->refine((i < 4 ? 4 : 10));
+ }
+ }
+
+ CouplingTimeInterface::solveInitialProblem(adaptInfo);
+ }
+
+protected:
+
+ CHNS_DrivenCavity *chnsProb;
+ PhaseFieldRefinement* refFunction;
+ RefinementLevelDOF *refinement;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+ AMDiS::init(argc, argv);
+
+ CHNS_DrivenCavity chnsProb("chns");
+ chnsProb.initialize(INIT_ALL);
+
+ RefinementTimeInterface timeInterface(&chnsProb);
+
+ // Adapt-Infos
+ AdaptInfo adaptInfo("adapt", chnsProb.getNumComponents());
+ ExtendedRosenbrockAdaptInstationary<CHNS_DrivenCavity> adaptInstat("adapt", chnsProb, adaptInfo, timeInterface, adaptInfo);
+
+
+ ptime start_time = microsec_clock::local_time();
+ chnsProb.initTimeInterface();
+ int error_code = adaptInstat.adapt();
+ time_duration td = microsec_clock::local_time()-start_time;
+
+ MSG("elapsed time= %d sec\n", td.total_seconds());
+
+ AMDiS::finalize();
+
+ return error_code;
+};
diff --git a/extensions/demo/src/fsi_explicit/ElasticityNavierStokes.h b/extensions/demo/src/fsi_explicit/ElasticityNavierStokes.h
index 9c120acb06878b58ebceda5e116ab5fc1254b12a..8a3e95ee0ee64d0391702bd5d41fb17797bfc3de 100644
--- a/extensions/demo/src/fsi_explicit/ElasticityNavierStokes.h
+++ b/extensions/demo/src/fsi_explicit/ElasticityNavierStokes.h
@@ -96,10 +96,10 @@ public:
addTimeInterface(elastProb2);
addTimeInterface(nsProb);
-// CouplingProblemStat::initialize(INIT_ALL | INIT_EXACT_SOLUTION);
- nsProb->initialize(INIT_ALL);
- elastProb->initialize(INIT_ALL);
- elastProb2->initialize(INIT_ALL);
+ CouplingProblemStat::initialize(INIT_ALL | INIT_EXACT_SOLUTION);
+// elastProb->initialize(INIT_ALL);
+// elastProb2->initialize(INIT_ALL);
+// nsProb->initialize(INIT_ALL - CREATE_MESH, elastProb2->getProblem(), INIT_MESH);
dim = elastProb->getMesh()->getDim();
initData();
@@ -223,18 +223,18 @@ public:
CouplingTimeInterface::closeTimestep(adaptInfo);
- flag->move(elastProb->getDeformation());
+// flag->move(elastProb->getDeformation());
updateMesh(elastProb->getDeformation(), parametricCoords1, parametric1);
updateMesh(elastProb2->getDeformation(), parametricCoords2, parametric2);
// refinement1->refine(5);
// refinement2->refine(5);
}
- void updateMesh(DOFVector<WorldVector<double> >* vec, WorldVector<DOFVector<double>*>& parametricCoords, ParametricFirstOrder* parametric)
+ void updateMesh(WorldVector<DOFVector<double>*> vec, WorldVector<DOFVector<double>*>& parametricCoords, ParametricFirstOrder* parametric)
{
FUNCNAME("ParametricSphere::buildAfterCoarsen()");
MSG("calculation of parametric coordinates\n");
- const FiniteElemSpace* feSpace = vec->getFeSpace();
+ const FiniteElemSpace* feSpace = vec[0]->getFeSpace();
int dim = feSpace->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
WorldVector<double> newCoords;
@@ -259,10 +259,9 @@ public:
basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
for (int i = 0; i < dim + 1; i++) {
dof = localIndices[i];
- newCoords = (*vec)[dof];
if (!visited[dof]) {
for (int j = 0; j < dow; j++)
- (*(parametricCoords[j]))[dof] += newCoords[j];
+ (*(parametricCoords[j]))[dof] += (*vec[j])[dof];
visited[dof] = true;
}
diff --git a/extensions/demo/src/fsi_explicit/fluidStructureInteraction.cc b/extensions/demo/src/fsi_explicit/fluidStructureInteraction.cc
index e8fba211a89baef3e742e1413d208ceb6803a7ef..c4476f575d4b60413b1ea1c6dbaa19a32794f639 100644
--- a/extensions/demo/src/fsi_explicit/fluidStructureInteraction.cc
+++ b/extensions/demo/src/fsi_explicit/fluidStructureInteraction.cc
@@ -22,10 +22,7 @@ struct NormalStress : TertiaryAbstractFunction<double, WorldVector<double>, Worl
double operator()(const WorldVector<double>& normal, const WorldVector<double>& grdU0, const WorldVector<double>& grdU1) const
{
- if (comp == 0)
- return viscosity * (2.0*grdU0[0]*normal[0] + (grdU1[0] + grdU0[1])*normal[1]);
- else
- return viscosity * (2.0*grdU1[1]*normal[1] + (grdU1[0] + grdU0[1])*normal[0]);
+ return viscosity * (2.0*grdU0[comp]*normal[comp] + (grdU1[0] + grdU0[1])*normal[1-comp]);
}
private:
int comp;
@@ -38,7 +35,7 @@ struct NormalStressP : TertiaryAbstractFunction<double, WorldVector<double>, dou
double operator()(const WorldVector<double>& normal, const double& stress, const double& p) const
{
- return stress-p*normal[comp];
+ return stress - p*normal[comp];
}
private:
int comp;
@@ -73,12 +70,12 @@ public:
typedef LinearElasticity super;
public:
- Elasticity_Obstacle(std::string name_) : super(name_), stress(NULL) {}
+ Elasticity_Obstacle(std::string name_) : super(name_) {}
~Elasticity_Obstacle()
{
- delete deformation;
- delete deformationVelocity;
+ for (size_t i = 0; i < dow; i++)
+ delete stress[i];
}
void initData()
@@ -86,19 +83,37 @@ public:
super::initData();
dim = getMesh()->getDim();
- stress = new DOFVector<WorldVector<double> >(getFeSpace(0), "stress");
- deformation = new DOFVector<WorldVector<double> >(getFeSpace(0), "deformation");
- deformationVelocity = new DOFVector<WorldVector<double> >(getFeSpace(dow), "deformation_velocity");
+ for (size_t i = 0; i < dow; i++) {
+ stress[i] = new DOFVector<double>(getFeSpace(0), "stress_"+Helpers::toString(i));
+ stress[i]->set(0.0);
+ }
+ for (size_t i = 0; i < dow; i++) {
+ deformation[i] = prob->getSolution()->getDOFVector(i);
+ deformationVelocity[i] = prob->getSolution()->getDOFVector(dow+i);
+ }
+
+ std::vector<DOFVector<double>*> vecs;
+ for (size_t i = 0; i < dow; i++)
+ vecs.push_back(stress[i]);
+ fileWriterStress = new FileWriter(name + "->stress->output", getFeSpace()->getMesh(), vecs);
}
- DOFVector<WorldVector<double> >* getDeformation()
+ DOFVector<double>* getDeformation(int i)
{
- return deformation;
+ return deformation[i];
+ }
+ WorldVector<DOFVector<double>*> getDeformation()
+ {
+ return deformation;
}
- DOFVector<WorldVector<double> >* getDeformationVelocity()
+ DOFVector<double>* getDeformationVelocity(int i)
+ {
+ return deformationVelocity[i];
+ }
+ WorldVector<DOFVector<double>*> getDeformationVelocity()
{
- return deformationVelocity;
+ return deformationVelocity;
}
void setNormalStress(DOFVector<WorldVector<double> >* stress_)
@@ -109,22 +124,20 @@ public:
void transferInitialSolution(AdaptInfo* adaptInfo)
{
super::transferInitialSolution(adaptInfo);
- calcDeformation();
- calcDeformationVelocity();
}
void initTimestep(AdaptInfo* adaptInfo)
{
super::initTimestep(adaptInfo);
TEST_EXIT(boundaryStress != NULL)("NULL-Pointer problem!\n");
- interpol_coords(*boundaryStress, *stress);
+ for (size_t i = 0; i < dow; i++)
+ transformDOF_coords(*boundaryStress, *stress[i], new Component(i));
+ fileWriterStress->writeFiles(adaptInfo, false);
}
void closeTimestep(AdaptInfo* adaptInfo)
{
super::closeTimestep(adaptInfo);
- calcDeformation();
- calcDeformationVelocity();
}
void fillBoundaryConditions()
@@ -135,42 +148,24 @@ public:
double zeroValue = 0.0;
// +--+
- // | | 10
- // 1| +==========
+ // | | 1
+ //10| +==========
// | |
// +--+
for (int i = 0; i < dow; i++) {
- prob->addDirichletBC(1, i, i, zero);
- prob->addNeumannBC(10, i, i, new VectorBC(i, stress));
+ prob->addDirichletBC(10, i, i, zero);
+ prob->addNeumannBC(1, i, i, stress[i]);
}
}
private:
- DOFVector<WorldVector<double> >* deformation;
- DOFVector<WorldVector<double> >* deformationVelocity;
- DOFVector<WorldVector<double> >* stress;
+ WorldVector<DOFVector<double>*> deformation;
+ WorldVector<DOFVector<double>*> deformationVelocity;
+ WorldVector<DOFVector<double>*> stress;
DOFVector<WorldVector<double> >* boundaryStress;
-
- void calcDeformation()
- {
- if (dow == 1)
- transformDOF(prob->getSolution()->getDOFVector(0), deformation, new AMDiS::Vec1WorldVec<double>);
- else if (dow == 2)
- transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), deformation, new AMDiS::Vec2WorldVec<double>);
- else if (dow == 3)
- transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), prob->getSolution()->getDOFVector(2), deformation, new AMDiS::Vec3WorldVec<double>);
- }
- void calcDeformationVelocity()
- {
- if (dow == 1)
- transformDOF(prob->getSolution()->getDOFVector(dow), deformationVelocity, new AMDiS::Vec1WorldVec<double>);
- else if (dow == 2)
- transformDOF(prob->getSolution()->getDOFVector(dow), prob->getSolution()->getDOFVector(dow+1), deformationVelocity, new AMDiS::Vec2WorldVec<double>);
- else if (dow == 3)
- transformDOF(prob->getSolution()->getDOFVector(dow), prob->getSolution()->getDOFVector(dow+1), prob->getSolution()->getDOFVector(dow+2), deformationVelocity, new AMDiS::Vec3WorldVec<double>);
- }
+ FileWriter* fileWriterStress;
};
///_______________________________________________________________________________________________________________________
@@ -181,11 +176,12 @@ public:
typedef StandardBaseProblem super;
public:
- Elasticity_Fluid(std::string name_) : super(name_), boundaryDeformation(NULL) {}
+ Elasticity_Fluid(std::string name_) : super(name_) {}
~Elasticity_Fluid()
{
- delete deformation;
+ for (size_t i = 0; i < dow; i++)
+ delete deformationVec[i];
}
void initData()
@@ -193,15 +189,24 @@ public:
super::initData();
dim = getMesh()->getDim();
- deformation = new DOFVector<WorldVector<double> >(getFeSpace(0), "deformation");
+ for (size_t i = 0; i < dow; i++) {
+ deformationVec[i] = new DOFVector<double>(getFeSpace(i), "deformation");
+ deformationVec[i]->set(0.0);
+ }
+ for (size_t i = 0; i < dow; i++)
+ deformation = prob->getSolution()->getDOFVector(i);
}
- DOFVector<WorldVector<double> >* getDeformation()
+ DOFVector<double>* getDeformation(int i)
+ {
+ return deformation[i];
+ }
+ WorldVector<DOFVector<double>*> getDeformation()
{
return deformation;
}
- void setDeformation(DOFVector<WorldVector<double> >* deformation_)
+ void setDeformation(WorldVector<DOFVector<double>*> deformation_)
{
boundaryDeformation = deformation_;
}
@@ -209,20 +214,19 @@ public:
void transferInitialSolution(AdaptInfo* adaptInfo)
{
super::transferInitialSolution(adaptInfo);
- calcDeformation();
}
void initTimestep(AdaptInfo* adaptInfo)
{
super::initTimestep(adaptInfo);
- TEST_EXIT(boundaryDeformation != NULL)("NULL-Pointer problem!\n");
- interpol_coords(*boundaryDeformation, *deformation);
+ TEST_EXIT(boundaryDeformation[0] != NULL)("NULL-Pointer problem!\n");
+ for (size_t i = 0; i < dow; i++)
+ interpol_coords(*boundaryDeformation[i], *deformationVec[i]);
}
void closeTimestep(AdaptInfo* adaptInfo)
{
super::closeTimestep(adaptInfo);
- calcDeformation();
}
void fillOperators()
@@ -243,34 +247,26 @@ public:
double zeroValue = 0.0;
// +--+
- // | | 10
- // 1| +==========
+ // | | 1
+ //10| +==========
// | |
// +--+
for (int i = 0; i < dow; i++) {
- prob->addDirichletBC(1, i, i, zero);
+ prob->addDirichletBC(10, i, i, zero);
prob->addDirichletBC(2, i, i, zero);
prob->addDirichletBC(3, i, i, zero);
prob->addDirichletBC(4, i, i, zero);
prob->addDirichletBC(5, i, i, zero);
- prob->addDirichletBC(10, i, i, new VectorBC(i, deformation));
+
+ prob->addDirichletBC(1, i, i, deformationVec[i]);
}
}
private:
- DOFVector<WorldVector<double> >* deformation;
- DOFVector<WorldVector<double> >* boundaryDeformation;
-
- void calcDeformation()
- {
- if (dow == 1)
- transformDOF(prob->getSolution()->getDOFVector(0), deformation, new AMDiS::Vec1WorldVec<double>);
- else if (dow == 2)
- transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), deformation, new AMDiS::Vec2WorldVec<double>);
- else if (dow == 3)
- transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), prob->getSolution()->getDOFVector(2), deformation, new AMDiS::Vec3WorldVec<double>);
- }
+ WorldVector<DOFVector<double>*> deformation; // solution
+ WorldVector<DOFVector<double>*> boundaryDeformation;
+ WorldVector<DOFVector<double>*> deformationVec; // boundaryDeformation->deformationVec
};
///_______________________________________________________________________________________________________________________
@@ -281,18 +277,30 @@ public:
typedef NavierStokes_TaylorHood super;
public:
- NS_Channel(std::string name_) : super(name_), boundaryVelocity(NULL) {}
+ NS_Channel(std::string name_) : super(name_) {}
~NS_Channel()
{
delete inflowBC;
+ delete fileWriterStress;
+ delete stress;
+ for (size_t i = 0; i < dow; i++)
+ delete velocityVec[i];
}
void initData()
{
super::initData();
stress = new DOFVector<WorldVector<double> >(getFeSpace(0), "stress");
+ WorldVector<double> zeroVec; zeroVec.set(0.0);
+ stress->set(zeroVec);
+ for (size_t i = 0; i < dow; i++) {
+ velocityVec[i] = new DOFVector<double>(getFeSpace(i), "velocity_i");
+ velocityVec[i]->set(0.0);
+ }
+
dim = getMesh()->getDim();
+ fileWriterStress = new FileVectorWriter(name + "->stress->output", getFeSpace()->getMesh(), stress);
}
void solveInitialProblem(AdaptInfo* adaptInfo)
@@ -301,13 +309,19 @@ public:
inflowBC->setTimePtr(adaptInfo->getTimePtr());
}
+
+ void transferInitialSolution(AdaptInfo* adaptInfo)
+ {
+ super::transferInitialSolution(adaptInfo);
+ fileWriterStress->writeFiles(adaptInfo, false);
+ }
DOFVector<WorldVector<double> >* getNormalStress()
{
return stress;
}
- void setDeformationVelocity(DOFVector<WorldVector<double> >* deformationVel_)
+ void setDeformationVelocity(WorldVector<DOFVector<double>*> deformationVel_)
{
boundaryVelocity = deformationVel_;
}
@@ -315,14 +329,16 @@ public:
void initTimestep(AdaptInfo* adaptInfo)
{
super::initTimestep(adaptInfo);
- TEST_EXIT(boundaryVelocity != NULL)("NULL-Pointer problem!\n");
- interpol_coords(*boundaryVelocity, *velocity);
+ TEST_EXIT(boundaryVelocity[0] != NULL)("NULL-Pointer problem!\n");
+ for (size_t i = 0; i < dow; i++)
+ interpol_coords(*boundaryVelocity[i], *velocityVec[i]);
}
void closeTimestep(AdaptInfo* adaptInfo)
{
super::closeTimestep(adaptInfo);
- calcNormalStress();
+ calcNormalStress(1);
+ fileWriterStress->writeFiles(adaptInfo, false);
}
void fillBoundaryConditions()
@@ -339,8 +355,8 @@ public:
/// at rigid wall: no-slip boundary condition
for (size_t i = 0; i < dow; i++) {
- getProblem(0)->addDirichletBC(1, i, i, zero);
- getProblem(0)->addDirichletBC(10, i, i, new VectorBC(i, velocity));
+ getProblem(0)->addDirichletBC(10, i, i, zero);
+ getProblem(0)->addDirichletBC(1, i, i, velocityVec[i]);
getProblem(0)->addDirichletBC(4, i, i, zero);
getProblem(0)->addDirichletBC(5, i, i, zero);
}
@@ -354,12 +370,14 @@ public:
/// at left wall: prescribed velocity
getProblem(0)->addDirichletBC(2, 0, 0, inflowBC);
getProblem(0)->addDirichletBC(2, 1, 1, zero);
+ /// at outflow boundary: pressure const 0.0
getProblem(0)->addDirichletBC(3, 1, 1, zero);
}
private:
DOFVector<WorldVector<double> >* stress;
- DOFVector<WorldVector<double> >* boundaryVelocity;
+ WorldVector<DOFVector<double>*> boundaryVelocity;
+ WorldVector<DOFVector<double>*> velocityVec; // boundaryVelocity->velocityVec
void calcNormalStress(BoundaryType boundaryType = 1)
{
@@ -374,35 +392,44 @@ private:
WorldVector<double> normal;
- std::map<DegreeOfFreedom, double> volume;
+ std::map<DegreeOfFreedom, bool> visited;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasisFcts = basFcts->getNumber();
- DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_DET | Mesh::FILL_COORDS);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
- mtl::dense_vector<double> localUh(nBasisFcts);
std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
+
+ // for pressure
+ const FiniteElemSpace* feSpaceP = getFeSpace(dow);
+ const BasisFunction *basFctsP = feSpaceP->getBasisFcts();
+ int nBasisFctsP = basFctsP->getNumber();
+ mtl::dense_vector<double> pLocalCoeffs(nBasisFctsP);
while (elInfo) {
for (int face = 0; face < dim + 1; face++) {
- if (elInfo->getBoundary(face) == boundaryType) {
+ if (elInfo->getBoundary(face) != 0) { //boundaryType) {
elInfo->getNormal(face, normal);
- double det = elInfo->getDet();
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+
+ // pressure
+ getSolution()->getDOFVector(dow)->getLocalVector(elInfo->getElement(), pLocalCoeffs);
for (int i = 0; i < nBasisFcts; i++) {
- double p = (*getSolution()->getDOFVector(dow))[localIndices[i]];
-
- for (int j = 0; j < Global::getGeo(WORLD); j++) {
- double tmp = NormalStress(j, viscosity)(normal, (*grdU[0])[localIndices[i]], (*grdU[1])[localIndices[i]]);
- double result = NormalStressP(j)(normal, tmp, p);
- (*stress)[localIndices[i]][j] += result * det;
+ if (!visited[localIndices[i]]) {
+
+ double p = basFctsP->evalUh(*basFcts->getCoords(i), pLocalCoeffs);
+
+ for (int j = 0; j < Global::getGeo(WORLD); j++) {
+ double tmp = NormalStress(j, density*viscosity)(normal, (*grdU[0])[localIndices[i]], (*grdU[1])[localIndices[i]]);
+ double result = NormalStressP(j)(normal, tmp, p);
+ (*stress)[localIndices[i]][j] = result;
+ }
+ visited[localIndices[i]] = true;
}
- volume[localIndices[i]] += det;
}
}
}
@@ -410,20 +437,13 @@ private:
elInfo = stack.traverseNext(elInfo);
}
- std::map<DegreeOfFreedom, double>::iterator volIt;
-
- for (volIt = volume.begin(); volIt != volume.end(); ++volIt) {
- for (int j = 0; j < Global::getGeo(WORLD); j++) {
- (*stress)[volIt->first][j] *= 1.0 / volIt->second;
- }
- }
-
for (size_t i = 0; i < dow; i++)
delete grdU[i];
}
// Polygon* obstacle;
InflowBC* inflowBC;
+ FileVectorWriter *fileWriterStress;
};
diff --git a/extensions/demo/src/movingMesh.cc b/extensions/demo/src/movingMesh.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a50e582a4a940b1d6e8662b55e09def95712073e
--- /dev/null
+++ b/extensions/demo/src/movingMesh.cc
@@ -0,0 +1,205 @@
+#include "AMDiS.h"
+#include "movingMesh.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+///_______________________________________________________________________________________________________________________
+
+class Elasticity_Fluid : public StandardBaseProblem
+{
+public:
+ typedef StandardBaseProblem super;
+
+public:
+ Elasticity_Fluid(std::string name_) : super(name_), displacement(NULL) {}
+
+ ~Elasticity_Fluid()
+ { }
+
+ void setDisplacement(AbstractFunction<double, WorldVector<double> >* displacement_)
+ {
+ displacement = displacement_;
+ }
+
+ void initData()
+ {
+ super::initData();
+ dim = getMesh()->getDim();
+
+ displacementDOF = new DOFVector<double>(getFeSpace(), "displacement_fluid");
+ }
+
+ void fillOperators()
+ {
+ // ===== create matrix operator =====
+ for (size_t i = 0; i < dow; i++) {
+ Operator *matrixOperator = new Operator(getFeSpace());
+ matrixOperator->addTerm(new Simple_SOT);
+ prob->addMatrixOperator(matrixOperator, i, i);
+
+// Operator *rhsOperator = new Operator(getFeSpace());
+// rhsOperator->addTerm(new CoordsAtQP_ZOT(new AMDiS::Const<double, WorldVector<double> >(1.0)));
+// prob->addVectorOperator(rhsOperator, i);
+ }
+ }
+
+ void fillBoundaryConditions()
+ {
+ MSG("Elasticity_Fluid::fillBoundaryConditions\n");
+ super::fillBoundaryConditions();
+
+ AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+ double zeroValue = 0.0;
+
+ // +--+
+ // | | 1
+ //10| +==========
+ // | |
+ // +--+
+
+ for (int i = 0; i < dow; i++) {
+ prob->addDirichletBC(10, i, i, zero);
+ prob->addDirichletBC(2, i, i, zero);
+ prob->addDirichletBC(3, i, i, zero);
+ prob->addDirichletBC(4, i, i, zero);
+ prob->addDirichletBC(5, i, i, zero);
+ }
+ prob->addDirichletBC(1, 0, 0, zero);
+ prob->addDirichletBC(1, 1, 1, displacementDOF);
+ }
+
+ void initTimestep(AdaptInfo *adaptInfo)
+ { MSG("initTimestep()\n");
+
+ super::initTimestep(adaptInfo);
+ displacementDOF->interpol(displacement);
+ }
+
+ void closeTimestep(AdaptInfo *adaptInfo)
+ { MSG("closeTimestep()\n");
+
+ super::closeTimestep(adaptInfo);
+ VtkWriter::writeFile(displacementDOF, "displacement_fluid.vtu");
+ }
+
+private:
+ AbstractFunction<double, WorldVector<double> >* displacement;
+ DOFVector<double>* displacementDOF;
+};
+
+
+///_______________________________________________________________________________________________________________________
+
+class Elasticity_Solid : public StandardBaseProblem
+{
+public:
+ typedef StandardBaseProblem super;
+
+public:
+ Elasticity_Solid(std::string name_) : super(name_), displacement(NULL) {}
+
+ ~Elasticity_Solid()
+ { }
+
+ void initData()
+ {
+ super::initData();
+ dim = getMesh()->getDim();
+
+ displacementDOF = new DOFVector<double>(getFeSpace(), "displacement_fluid");
+ }
+
+ void setDisplacement(AbstractFunction<double, WorldVector<double> >* displacement_)
+ {
+ displacement = displacement_;
+ }
+
+ void fillOperators()
+ {
+ // ===== create matrix operator =====
+ for (size_t i = 0; i < dow; i++) {
+ Operator *matrixOperator = new Operator(getFeSpace());
+ matrixOperator->addTerm(new Simple_SOT);
+ prob->addMatrixOperator(matrixOperator, i, i);
+
+// Operator *rhsOperator = new Operator(getFeSpace());
+// rhsOperator->addTerm(new CoordsAtQP_ZOT(new AMDiS::Const<double, WorldVector<double> >(1.0)));
+// prob->addVectorOperator(rhsOperator, i);
+ }
+ }
+
+ void fillBoundaryConditions()
+ {
+ MSG("Elasticity_Solid::fillBoundaryConditions\n");
+ super::fillBoundaryConditions();
+
+ AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+ double zeroValue = 0.0;
+
+ // +--+
+ // | | 1
+ //10| +==========
+ // | |
+ // +--+
+
+ for (int i = 0; i < dow; i++) {
+ prob->addDirichletBC(10, i, i, zero);
+ prob->addDirichletBC(2, i, i, zero);
+ prob->addDirichletBC(3, i, i, zero);
+ prob->addDirichletBC(4, i, i, zero);
+ prob->addDirichletBC(5, i, i, zero);
+ }
+ prob->addDirichletBC(1, 0, 0, zero);
+ prob->addDirichletBC(1, 1, 1, displacementDOF);
+ }
+
+ void initTimestep(AdaptInfo *adaptInfo)
+ { MSG("initTimestep()\n");
+
+ super::initTimestep(adaptInfo);
+ displacementDOF->interpol(displacement);
+ }
+
+ void closeTimestep(AdaptInfo *adaptInfo)
+ { MSG("closeTimestep()\n");
+
+ super::closeTimestep(adaptInfo);
+ VtkWriter::writeFile(displacementDOF, "displacement_solid.vtu");
+ }
+
+private:
+ AbstractFunction<double, WorldVector<double> >* displacement;
+ DOFVector<double>* displacementDOF;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+ AMDiS::init(argc, argv);
+
+ Elasticity_Fluid elastProb("elasticity->fluid");
+ Elasticity_Solid elastProb2("elasticity->solid");
+
+ MovingMesh<Elasticity_Fluid, Elasticity_Solid> mainProb("main", &elastProb, &elastProb2);
+
+ // Adapt-Infos
+ AdaptInfo adaptInfo("adapt", mainProb.getNumComponents());
+ mainProb.initialize(&adaptInfo);
+
+ // adaption loop - solve ch-prob and ns-prob
+ AdaptInstationary adaptInstat("adapt", mainProb, adaptInfo, mainProb, adaptInfo);
+
+ ptime start_time = microsec_clock::local_time();
+ int error_code = adaptInstat.adapt();
+ time_duration td = microsec_clock::local_time()-start_time;
+
+ MSG("elapsed time= %d sec\n", td.total_seconds());
+
+ AMDiS::finalize();
+
+ return error_code;
+};
diff --git a/extensions/demo/src/movingMesh.h b/extensions/demo/src/movingMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ab92d470af1467fec3ca018f549dc31b605ad0c
--- /dev/null
+++ b/extensions/demo/src/movingMesh.h
@@ -0,0 +1,376 @@
+/** \file MovingMesh.h */
+
+#ifndef MOVING_MESH_H
+#define MOVING_MESH_H
+
+#include "BaseProblem.h"
+
+// coupling structures
+#include "CouplingIterationInterface.h"
+#include "CouplingTimeInterface.h"
+#include "CouplingProblemStat.h"
+
+// structures for local refinement
+// #include "Refinement.h"
+// #include "MeshFunction_Level.h"
+
+// #include "SignedDistFunctors.h"
+// #include "PhaseFieldConvert.h"
+
+using namespace AMDiS;
+
+
+struct MinWrapper : AbstractFunction<double, WorldVector<double> >
+{
+ MinWrapper(AbstractFunction<double, WorldVector<double> >* dist1_, AbstractFunction<double, WorldVector<double> >* dist2_)
+ : dist1(dist1_), dist2(dist2_) {}
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return std::min((*dist1)(x), (*dist2)(x));
+ }
+private:
+ AbstractFunction<double, WorldVector<double> >* dist1;
+ AbstractFunction<double, WorldVector<double> >* dist2;
+};
+
+
+struct BeamDisplacement1d : AbstractFunction<double, double>
+{
+ BeamDisplacement1d(double P_, double L_, double EI_, double* time_ = NULL, double* tau_ = NULL) : P(P_), L(L_), EI(EI_), time(time_), tau(tau_) {}
+
+ double operator()(const double& x) const
+ {
+ return (*time < DBL_TOL ? 0.0 : (sin(*time)-sin(*time - *tau))*P*sqr(x)*(3.0*L-x)/(6.0*EI));
+ }
+
+ void setTimePtr(double* time_)
+ {
+ time = time_;
+ }
+
+ void setTauPtr(double* tau_)
+ {
+ tau = tau_;
+ }
+
+private:
+ double P;
+ double L;
+ double EI;
+ double* time;
+ double* tau;
+};
+
+struct BeamDisplacement : AbstractFunction<double, WorldVector<double> >
+{
+ BeamDisplacement(double P, double L, double EI, WorldVector<double> fixedPoint_, double* time = NULL, double* tau = NULL)
+ : fct(new BeamDisplacement1d(P, L, EI, time, tau)), fixedPoint(fixedPoint_) {}
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ WorldVector<double> displacement; displacement[0] = 0.0;
+ displacement[1] = (x[0] >= fixedPoint[0] ? (*fct)(x[0] - fixedPoint[0]) : 0.0);
+ return displacement[1];
+ }
+
+ void setTimePtr(double* time_)
+ {
+ fct->setTimePtr(time_);
+ }
+
+ void setTauPtr(double* tau_)
+ {
+ fct->setTauPtr(tau_);
+ }
+
+private:
+ BeamDisplacement1d* fct;
+ WorldVector<double> fixedPoint;
+};
+
+/**
+ * \ingroup Problem
+ *
+ * \brief
+ */
+template<typename Elasticity_Type1 = StandardBaseProblem,
+ typename Elasticity_Type2 = StandardBaseProblem>
+class MovingMesh : public CouplingIterationInterface,
+ public CouplingTimeInterface,
+ public CouplingProblemStat
+{
+public:
+
+ MovingMesh(std::string name_, Elasticity_Type1 *elastProb_, Elasticity_Type2 *elastProb2_)
+ : CouplingProblemStat(name_),
+ elastProb(elastProb_),
+ elastProb2(elastProb2_)
+// refFunction1(NULL),
+// refFunction2(NULL),
+// refinement1(NULL),
+// refinement2(NULL)
+ {
+ dow = Global::getGeo(WORLD);
+ }
+
+
+ ~MovingMesh() {
+// if (refFunction1 != NULL)
+// delete refFunction1;
+// if (refinement1 != NULL)
+// delete refinement1;
+// if (refFunction2 != NULL)
+// delete refFunction2;
+// if (refinement2 != NULL)
+// delete refinement2;
+
+ for (int i = 0; i < Global::getGeo(WORLD); i++) {
+ delete parametricCoords1[i];
+ delete parametricCoords2[i];
+ }
+
+ delete parametric1;
+ delete parametric2;
+ }
+
+
+ /**
+ * Add the problems to the iterationInterface, timeInterface and couplingProblemStat.
+ * As a consequence all problem can be initialized one after another and in the
+ * adaption loop they are solved in rotation.
+ * At the end the problems are filled with operators and coupling operators as well as
+ * boundary conditions are added.
+ *
+ * In the adaption loop the problems are solved the same order as they are added to the
+ * iterationInterface in this method. This order can be changed manually in the oneIteration
+ * method.
+ **/
+ void initialize(AdaptInfo *adaptInfo)
+ {
+ for (size_t i = 0; i < elastProb->getNumProblems(); i++)
+ addProblem(elastProb->getProblem(i));
+ for (size_t i = 0; i < elastProb2->getNumProblems(); i++)
+ addProblem(elastProb2->getProblem(i));
+
+ addIterationInterface(elastProb);
+ addIterationInterface(elastProb2);
+
+ addTimeInterface(elastProb);
+ addTimeInterface(elastProb2);
+
+// CouplingProblemStat::initialize(INIT_ALL | INIT_EXACT_SOLUTION);
+ elastProb->initialize(INIT_ALL);
+ elastProb2->initialize(INIT_ALL);
+ dim = elastProb->getMesh()->getDim();
+
+ initData();
+ elastProb->initData();
+ elastProb2->initData();
+ fillCouplingBoundaryConditions();
+
+ elastProb->fillOperators();
+ elastProb2->fillOperators();
+
+ elastProb->fillBoundaryConditions();
+ elastProb2->fillBoundaryConditions();
+ }
+
+
+ void initData()
+ { MSG("initData()\n");
+
+ // ===== create coords vector =====
+ for (int i = 0; i < Global::getGeo(WORLD); i++) {
+ parametricCoords1[i] = new DOFVector<double>(elastProb->getFeSpace(0),
+ "parametric1 coords "+boost::lexical_cast<std::string>(i));
+ parametricCoords2[i] = new DOFVector<double>(elastProb2->getFeSpace(0),
+ "parametric2 coords "+boost::lexical_cast<std::string>(i));
+ }
+
+ // ===== create parametric object =====
+ parametric1 = new ParametricFirstOrder(¶metricCoords1);
+ parametric2 = new ParametricFirstOrder(¶metricCoords2);
+
+ double P = 0.2,L = 0.46,EI = 1.0;
+ Parameters::get("beam->P", P);
+ Parameters::get("beam->L", L);
+ Parameters::get("beam->EI", EI);
+ WorldVector<double> fixedPoint;
+ Parameters::get("beam->fixedPoint", fixedPoint);
+ beamDisplacement = new BeamDisplacement(P, L, EI, fixedPoint);
+ }
+
+
+ void fillCouplingBoundaryConditions()
+ { MSG("fillCouplingBoundaryConditions()\n");
+ elastProb->setDisplacement(beamDisplacement);
+ elastProb2->setDisplacement(beamDisplacement);
+ }
+
+
+ /**
+ * Solves the initial problems.Before this is done the mesh is refined by the Refinement-class
+ * This provides local refinement functions for phaseField-refinement or signed-distance-refinement.
+ * See MeshFunction_Level.h for some details.
+ **/
+ void solveInitialProblem(AdaptInfo *adaptInfo)
+ { MSG("solveInitialProblem()\n");
+
+ beamDisplacement->setTimePtr(adaptInfo->getTimePtr());
+ beamDisplacement->setTauPtr(adaptInfo->getTimestepPtr());
+
+// size_t nVertices = 0;
+// WorldVector<double> c,c2;
+// Parameters::get("obstacle->num vertices",nVertices);
+// std::vector<WorldVector<double> > v(nVertices,c);
+// for (size_t i = 0; i < nVertices; i++)
+// Parameters::get("obstacle->vertex["+boost::lexical_cast<std::string>(i)+"]",v[i]);
+// v.push_back(v[0]);
+//
+// obstacle = new Polygon(v);
+//
+// Parameters::get("flag->num vertices",nVertices);
+// std::vector<WorldVector<double> > v2(nVertices,c);
+// for (size_t i = 0; i < nVertices; i++)
+// Parameters::get("flag->vertex["+boost::lexical_cast<std::string>(i)+"]",v2[i]);
+// v2.push_back(v2[0]);
+//
+// flag = new Polygon(v2);
+//
+// Parameters::get("connection->num vertices",nVertices);
+// std::vector<WorldVector<double> > v3(nVertices,c);
+// for (size_t i = 0; i < nVertices; i++)
+// Parameters::get("connection->vertex["+boost::lexical_cast<std::string>(i)+"]",v3[i]);
+// v3.push_back(v3[0]);
+//
+// signedDistFct = new MinWrapper(new Polygon(v2), new MinWrapper(obstacle,flag));
+//
+// refFunction1 = new SignedDistRefinement(elastProb->getMesh());
+// refinement1 = new RefinementLevelCoords2(
+// elastProb->getFeSpace(0),
+// refFunction1,
+// signedDistFct);
+//
+// refFunction2 = new SignedDistRefinement(elastProb2->getMesh());
+// refinement2 = new RefinementLevelCoords2(
+// elastProb2->getFeSpace(0),
+// refFunction2,
+// signedDistFct);
+
+ // initial refinement
+// refinement1->refine(10);
+// refinement2->refine(10);
+// flag->refine(10);
+
+ // solve all initial problems of the problems added to the CouplingTimeInterface
+ CouplingTimeInterface::solveInitialProblem(adaptInfo);
+
+ DOFVector<WorldVector<double> > coords1(parametricCoords1[0]->getFeSpace(), "coords");
+ parametricCoords1[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords1[0]->getFeSpace(), coords1);
+ DOFVector<WorldVector<double> > coords2(parametricCoords2[0]->getFeSpace(), "coords");
+ parametricCoords2[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords2[0]->getFeSpace(), coords2);
+ for (int i = 0; i < Global::getGeo(WORLD); i++) {
+ transformDOF(&coords1, parametricCoords1[i], new Component(i));
+ transformDOF(&coords2, parametricCoords2[i], new Component(i));
+ }
+
+ // ===== enable parametric traverse =====
+ elastProb->getMesh()->setParametric(parametric1);
+ elastProb2->getMesh()->setParametric(parametric2);
+ }
+
+
+ /**
+ * Called at the end of each timestep.
+ * If the phase-field has changed the mesh is updated by the refinement-method
+ **/
+ void closeTimestep(AdaptInfo *adaptInfo)
+ { MSG("closeTimestep()\n");
+
+ CouplingTimeInterface::closeTimestep(adaptInfo);
+
+ DOFVector<double> displacement(elastProb->getFeSpace(), "displacement");
+ displacement.interpol(beamDisplacement);
+ VtkWriter::writeFile(&displacement, "displacement.vtu");
+
+// flag->move(elastProb->getDeformation());
+ updateMesh(elastProb->getSolution(), parametricCoords1, parametric1);
+ updateMesh(elastProb2->getSolution(), parametricCoords2, parametric2);
+// refinement1->refine(5);
+// refinement2->refine(5);
+ }
+
+ void updateMesh(SystemVector* vec, WorldVector<DOFVector<double>*>& parametricCoords, ParametricFirstOrder* parametric)
+ {
+ FUNCNAME("ParametricSphere::buildAfterCoarsen()");
+ MSG("calculation of parametric coordinates\n");
+ const FiniteElemSpace* feSpace = vec->getFeSpace(0);
+ int dim = feSpace->getMesh()->getDim();
+ int dow = Global::getGeo(WORLD);
+ WorldVector<double> newCoords;
+
+ DegreeOfFreedom dof;
+ WorldVector<double> x;
+
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+ int numBasFcts = basFcts->getNumber();
+ std::vector<DegreeOfFreedom> localIndices(numBasFcts);
+ DOFAdmin *admin = feSpace->getAdmin();
+
+ // ===== disable parametric traverse =====
+ feSpace->getMesh()->setParametric(NULL);
+
+ std::map<DegreeOfFreedom, bool> visited;
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ while (elInfo) {
+ basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
+ for (int i = 0; i < dim + 1; i++) {
+ dof = localIndices[i];
+ if (!visited[dof]) {
+ for (int j = 0; j < dow; j++)
+ (*(parametricCoords[j]))[dof] += (*vec->getDOFVector(j))[dof];
+
+ visited[dof] = true;
+ }
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ // ===== enable parametric traverse =====
+ feSpace->getMesh()->setParametric(parametric);
+ }
+
+protected:
+
+ Elasticity_Type1 *elastProb; // LinearElasticity baseProblem
+ Elasticity_Type2 *elastProb2; // LinearElasticity baseProblem
+
+ /// DOFVector storing parametric coordinates.
+ WorldVector<DOFVector<double>*> parametricCoords1;
+ WorldVector<DOFVector<double>*> parametricCoords2;
+
+ /// Parametrizes vertex coordinates while mesh traversal.
+ ParametricFirstOrder *parametric1;
+ ParametricFirstOrder *parametric2;
+
+// Polygon* flag;
+// Polygon* obstacle;
+// AbstractFunction<double, WorldVector<double> >* signedDistFct;
+
+// SignedDistRefinement* refFunction1;
+// SignedDistRefinement* refFunction2;
+// RefinementLevelCoords2* refinement1;
+// RefinementLevelCoords2* refinement2;
+
+ BeamDisplacement* beamDisplacement;
+
+ unsigned dim; // dimension of the mesh
+ unsigned dow; // dimension of the world
+};
+
+
+#endif // MOVING_MESH_H
diff --git a/extensions/demo/src/navierStokes.h b/extensions/demo/src/navierStokes.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e1750459d7c2a2822716e92290b10b77e3c3e40
--- /dev/null
+++ b/extensions/demo/src/navierStokes.h
@@ -0,0 +1,43 @@
+/** \file navierStokes.h */
+
+#ifndef NAVIER_STOKES_H
+#define NAVIER_STOKES_H
+
+#include "AMDiS.h"
+#include "GeometryTools.h"
+
+struct InflowBC : AbstractFunction<double, WorldVector<double> >
+{
+ InflowBC(double H_=4.1, double Um_=1.5) : H(H_), Um(Um_) {}
+ double operator()(const WorldVector<double> &x) const {
+ double vel = 4.0 * Um * x[1] * (H - x[1]) / sqr(H);
+ return (*time < 2.0 ? vel*(1.0 - cos(m_pi * (*time)/2.0))/2.0 : vel);
+ }
+
+ void setTimePtr(double* time_)
+ {
+ time = time_;
+ }
+
+protected:
+ double H;
+ double Um;
+ double* time;
+};
+
+
+struct MinWrapper : AbstractFunction<double, WorldVector<double> >
+{
+ MinWrapper(AbstractFunction<double, WorldVector<double> >* dist1_, AbstractFunction<double, WorldVector<double> >* dist2_)
+ : dist1(dist1_), dist2(dist2_) {}
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return std::min((*dist1)(x), (*dist2)(x));
+ }
+private:
+ AbstractFunction<double, WorldVector<double> >* dist1;
+ AbstractFunction<double, WorldVector<double> >* dist2;
+};
+
+#endif