Commit 13a42661 authored by Praetorius, Simon's avatar Praetorius, Simon

extension

parent f042083c
......@@ -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,
......
#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;
};
#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;
};
......@@ -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;
}
......
......@@ -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];