Skip to content
Snippets Groups Projects
Commit 3064059b authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

demo directory cleaned up

parent bc2410a6
No related branches found
No related tags found
No related merge requests found
/** \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 {
return 4.0 * Um * x[1] * (H - x[1]) / sqr(H);
}
protected:
double H;
double Um;
};
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
#include "AMDiS.h"
#include "NavierStokesPhase_TaylorHood.h"
#include "navierStokes.h"
// #include "time/ExtendedRosenbrockAdaptInstationary.h"
#include "Refinement.h"
#include "MeshFunction_Level.h"
#include "PhaseFieldConvert.h"
#include "boost/date_time/posix_time/posix_time.hpp"
using namespace AMDiS;
using namespace boost::posix_time;
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<WorldVector<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_) {}
WorldVector<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;
}
void setTimePtr(double* time_)
{
fct->setTimePtr(time_);
}
void setTauPtr(double* tau_)
{
fct->setTauPtr(tau_);
}
private:
BeamDisplacement1d* fct;
WorldVector<double> fixedPoint;
};
class NS_Channel : public NavierStokesPhase_TaylorHood
{
public:
typedef NavierStokesPhase_TaylorHood super;
public:
NS_Channel(std::string name_) : super(name_) {}
~NS_Channel()
{
delete phaseField;
}
void initData()
{ FUNCNAME("NS_Channel::initData()");
super::initData();
phaseField = new DOFVector<double>(getFeSpace(0), "phaseField");
phaseField->set(1.0);
super::setPhase(phaseField);
}
void solveInitialProblem(AdaptInfo *adaptInfo)
{ FUNCNAME("NS_Channel::solveInitialProblem()");
super::solveInitialProblem(adaptInfo);
size_t nVertices = 0;
WorldVector<double> x;
Parameters::get("obstacle->num vertices",nVertices);
std::vector<WorldVector<double> > v(nVertices,x);
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);
refFunction = new SignedDistRefinement(getMesh());
refinement = new RefinementLevelCoords2(
getFeSpace(),
refFunction,
obstacle);
// initial refinement
refinement->refine(10);
phaseField->interpol(new SignedDistFctToPhaseField(getEpsilon(), obstacle, -3.0));
beamDisplacement->setTimePtr(adaptInfo->getTimePtr());
beamDisplacement->setTauPtr(adaptInfo->getTimestepPtr());
obstacle->refine(10);
}
protected:
void fillBoundaryConditions()
{ FUNCNAME("NS_Channel::fillBoundaryConditions()");
AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
size_t dow = Global::getGeo(WORLD);
double P = 0.1, L = 2.0, EI = 1.0;
WorldVector<double> fixedPoint; fixedPoint[0] = 2.5; fixedPoint[1] = 2.0;
Parameters::get("obstacle->P",P);
Parameters::get("obstacle->L",L);
Parameters::get("obstacle->EI",EI);
Parameters::get("obstacle->fixed point",fixedPoint);
beamDisplacement = new BeamDisplacement(P, L, EI, fixedPoint);
super::setBcFct(beamDisplacement);
super::fillBoundaryConditions();
// +------ 5 ------+
// | |
// 2 # <--1 3
// | |
// +------ 4 ------+
/// at rigid wall: no-slip boundary condition
for (size_t i = 0; i < dow; i++) {
getProblem(0)->addDirichletBC(4, i, i, zero);
getProblem(0)->addDirichletBC(5, i, i, zero);
}
/// dirichlet bc for pressure at one DOF
// getProblem(0)->addSingularDirichletBC(0, dow, dow, *zero);
double H = 4.1;
double Um = 1.5;
Parameters::get("mesh->H",H);
Parameters::get("ns->Um",Um);
/// at left wall: prescribed velocity
// getProblem(0)->addDirichletBC(2, 0, 0, new InflowBC(H,Um));
// getProblem(0)->addDirichletBC(2, 1, 1, zero);
}
void initTimestep(AdaptInfo* adaptInfo)
{
super::initTimestep(adaptInfo);
obstacle->move(beamDisplacement);
refinement->refine(5);
phaseField->interpol(new SignedDistFctToPhaseField(getEpsilon(), obstacle, -3.0));
}
private:
DOFVector<double>* phaseField;
BeamDisplacement* beamDisplacement;
Polygon* obstacle;
SignedDistRefinement* refFunction;
RefinementLevelCoords2* refinement;
};
int main(int argc, char** argv)
{ FUNCNAME("main");
AMDiS::init(argc, argv);
NS_Channel nsProb("ns");
nsProb.initialize(INIT_ALL | INIT_EXACT_SOLUTION);
// Adapt-Infos
AdaptInfo adaptInfo("adapt", nsProb.getNumComponents());
// adaption loop - solve ch-prob and ns-prob
AdaptInstationary adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
// ExtendedRosenbrockAdaptInstationary<NS_Channel> adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
// scale Mesh
WorldVector<double> scale; scale[0] = 25.0; scale[1] = 4.1;
Helpers::scaleMesh(nsProb.getMesh(), scale);
ptime start_time = microsec_clock::local_time();
nsProb.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;
};
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment