Skip to content
Snippets Groups Projects
Commit 4009c940 authored by Lars von Wolff's avatar Lars von Wolff
Browse files

removed archive

parent f1155d79
No related branches found
No related tags found
No related merge requests found
#pragma once
//#include "../myincludes.hh"
#include <cmath>
template<typename GV, typename RF, std::size_t dim_range>
class ZeroVectorFunction :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,
ZeroVectorFunction<GV,RF,dim_range> >,
public Dune::PDELab::InstationaryFunctionDefaults
{
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroVectorFunction> BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
ZeroVectorFunction(const GV & gv) : BaseT(gv) {}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y=0;
}
};
template<typename GV, typename RF>
class ZeroScalarFunction
: public ZeroVectorFunction<GV,RF,1>
{
public:
ZeroScalarFunction(const GV & gv) : ZeroVectorFunction<GV,RF,1>(gv) {}
};
template<typename GV, typename RF, int dim>
class PhiZero :
public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, PhiZero<GV,RF,dim> >
{
const RF mean;
const RF maxnoise;
const RF eps;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiZero<GV,RF,dim> > BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
PhiZero(const GV & gv, const RF& mean_, const RF& maxnoise_, const RF& eps_) :
BaseT(gv), mean(mean_), maxnoise(maxnoise_), eps(eps_) {}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
//Random Noise
//y[0] = mean+ maxnoise*(-1 + 2*((RF)rand())/RAND_MAX);
//Deterministic "pseudeorandom"
//y[0] = mean+ maxnoise* sin(sqrt(3)*100000*x[0]+sqrt(2))*cos(sqrt(2)*1000000*x[1]+sqrt(5))
//Bubble
//RF radius = sqrt( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5) );
//y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*(radius-0.2)));
//Stripes
// RF d = 0;
// if(x[0]<3./8)
// d = x[0]-1./4;
// else if(x[0] < 5./8)
// d = 1./8 - (x[0]-3./8);
// else
// d = -1./8 + (x[0]-5./8);
// y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
//Left/Right
//y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-0.75)));
//Top/Bottom (Rayleigh-Taylor Instability)
y[0] = 0.5*(1-tanh(1/eps/sqrt(2)*(x[1]-0.75-0.01*cos(x[0]*2*3.1415))));
}
};
template<typename GV, typename RF, int dim>
class ThreephasePhiZero :
public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, ThreephasePhiZero<GV,RF,dim> >
{
const RF mean;
const RF maxnoise;
const RF eps;
const int phaseNumber;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,ThreephasePhiZero<GV,RF,dim> > BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
ThreephasePhiZero(const GV & gv, const int& num_, const RF& mean_, const RF& maxnoise_, const RF& eps_) :
BaseT(gv), mean(mean_), maxnoise(maxnoise_), eps(eps_), phaseNumber(num_) {}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
//Stripes
// if(phaseNumber == 0)
// {
// RF d = (x[0]-3./4);
// y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
// }
// if(phaseNumber == 1)
// {
// RF d = -(x[0]-1./4);
// y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
// }
//Contact point
RF d1 = x[0]-0.5;
RF yx = 0.5*(1+tanh(1/eps/sqrt(2)*d1));
RF d2 = x[1]-0.75;
if(phaseNumber == 0)
d2 = -d2;
y[0]=0.5*(1+tanh(1/eps/sqrt(2)*d2))*yx;
}
};
template<typename GV, typename RF, int dim>
class UZero :
public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, UZero<GV,RF,dim> >
{
const RF mean;
const RF maxnoise;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,UZero<GV,RF,dim> > BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
UZero(const GV & gv, const RF& mean_, const RF& maxnoise_) :
BaseT(gv), mean(mean_), maxnoise(maxnoise_) {}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y[0] = mean;
//if (x[0]> 1-1e-5)
//{
// y[0]=mean-0.3;
//}
}
};
class U_BC : public Dune::PDELab::DirichletConstraintsParameters
{
public:
template<typename I>
bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
{
Dune::FieldVector<typename I::ctype, I::coorddimension>
xg = intersection.geometry().global( coord );
if (xg[0]> 1-1e-5)
{
return true;
}
return false;
}
};
class SinglePointPressureBC : public Dune::PDELab::DirichletConstraintsParameters
{
public:
template<typename I> bool isDirichlet(const I & intersection,
const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
{
Dune::FieldVector<typename I::ctype, I::coorddimension>
xg = intersection.geometry().global( coord );
//std::cout<< "x0: " << xg[0] << " x1: " << xg[1] << std::endl;
if (xg[0]< 1e-2)
{
//std::cout<< "Found Pressure Zero Point";
//return true;
}
return false;
}
};
template<typename GV, typename RF, int dim>
class IniVelocity :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
IniVelocity<GV,RF,dim> >
{
const RF meanflow_;
const RF timemax;
RF time;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,IniVelocity<GV,RF,dim> > BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
IniVelocity(const GV & gv, const RF& meanflow, const RF& tm) :
BaseT(gv)
, meanflow_(meanflow), timemax(tm)
{ time=0.0;}
template <typename T>
void setTime(T t){
time = t;
}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
RF radius = sqrt( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5) );
RF r = meanflow_*(0.5-radius);
if(r< 0)
{
r=0;
}
//y[0] = -r*(x[1]-0.5) * 1e-4; //Rotating flow
//y[1] = r*(x[0]-0.5) * 1e-4;
//y[1]=x[0]; //Shear flow
//y[0]=0;
//y[0]=0; //No flow
//y[1]=0;
RF tscale=fmin(time/timemax,1);
y[0]=0; //Shear flow boundary
if(x[0]>0.9999)
{
//y[1]=4*meanflow_*x[1]*(1-x[1])*tscale;
if (x[1]>0 && x[1]<0.75)
{
y[1]=64/9*meanflow_*x[1]*(0.75-x[1])*tscale;
}
}
}
};
template<typename GV, typename RF, int dim>
class TUPressure :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
TUPressure<GV,RF,dim> >
{
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TUPressure<GV,RF,dim> > BaseT;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
TUPressure(const GV & gv) :
BaseT(gv)
{ }
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
//y[0]=0.1*x[0];
}
template <typename T> void setTime(T t){}
};
#pragma once
//#include "../myincludes.hh"
#include <cmath>
template<typename Number>
class Problem
{
Number t;
const Number TurningPointLeft = 0.5 - 1/sqrt(12);
const Number TurningPointRight = 0.5 + 1/sqrt(12);
public:
typedef Number value_type;
Number eps, delta;
Number D;
Number dt;
Number uAst;
Number u_regularization;
Number mu;
Number rho_0;
Number rho_1;
Number rho_2;
Number M;
Number M_reduced;
Number fscaling;
Number SurfaceTensionScaling;
//Number rho (const Number& x) {return rho_0;}
Number rho (const Number& x){return x*rho_2 + (1-x)*rho_1;}
Number Sigma1,Sigma2,Sigma3,SigmaT;
// Constructor takes parameter
Problem (const Number& dt_, const Dune::ParameterTree& ptree)
: t(0.0)
{
eps = ptree.get("problem.eps",(Number)0.1);
delta = ptree.get("problem.delta",(Number)0.1);
M = ptree.get("problem.Mobility",(Number)0.1);
M_reduced = ptree.get("problem.MobilityReduced",(Number)0.1);
u_regularization = ptree.get("problem.URegularization",(Number)0.01);
D = ptree.get("physics.D",(Number)0.1);
uAst = ptree.get("physics.uAst",(Number)0.1);
rho_0 = ptree.get("physics.rho",(Number)1);
mu = ptree.get("physics.mu",(Number)0.001);
fscaling = ptree.get("phsyics.fscaling",(Number)1);
SurfaceTensionScaling = ptree.get("physics.SurfaceTensionScaling",(Number)1);
rho_1 = ptree.get("3p.rho1",(Number)1);
rho_2 = ptree.get("3p.rho2",(Number)1);
Sigma1 = ptree.get("3p.Sigma1",(Number)1);
Sigma2 = ptree.get("3p.Sigma2",(Number)1);
Sigma3 = ptree.get("3p.Sigma3",(Number)1);
SigmaT = 1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3);
dt=dt_;
}
Number Mob (const Number& u) const
{
return M; //M_reduced + (M-M_reduced)*q(u);
}
//Double Well
Number W (const Number& u) const
{
return u*u*(1-u)*(1-u);
}
// Double Well Derivative(/Secant)
Number dW (const Number& x, const Number& y) const
{
return (x+y)*(1 +x*x +y*y) -2*(x*x +x*y +y*y);
}
// Double Well Derivative
Number dW (const Number& x) const
{
return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
}
// Double Well Derivative of Convex Part
Number dWConvex (const Number& x) const
{
if(x<TurningPointLeft)
return dW(x)-dW(TurningPointLeft);
if(x>TurningPointRight)
return dW(x)-dW(TurningPointRight);
return 0;
}
// Double Well Derivative of Convcarve Part
Number dWConcarve (const Number& x) const
{
if(x<TurningPointLeft)
return dW(TurningPointLeft);
if(x>TurningPointRight)
return dW(TurningPointRight);
return dW(x);
}
Number dWdpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma1-SigmaT)*dW(pf1)-SigmaT*dW(pf2)-SigmaT*dW(pf3);
}
Number dWdpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma2-SigmaT)*dW(pf2)-SigmaT*dW(pf1)-SigmaT*dW(pf3);
}
Number dWdpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma3-SigmaT)*dW(pf3)-SigmaT*dW(pf1)-SigmaT*dW(pf2);
}
Number f (const Number& u) const
{
if(u<0)
return fscaling*(-1.0+1000*u);
else
return fscaling*(u*u-1);
}
Number q (const Number& pf) const
{
if (pf<0 || pf>1)
return 0;
return pf*(1-pf);
}
Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
{
if (pf1<0 || pf3<0)
return 0;
return pf1*pf3;
}
Number Diffusion(const Number& pf_c)
{
return D*fmax(1e-3,(pf_c-0.05)/0.95);
}
Number VelDissipation(const Number& pf_f)
{
return 10*(1-pf_f)*(1-pf_f);
}
//! Set time in instationary case
void setTime (Number t_)
{
t = t_;
}
};
#pragma once
//#include "../myincludes.hh"
#include <cmath>
template<typename Number>
class Problem
{
Number t;
const Number TurningPointLeft = 0.5 - 1/sqrt(12);
const Number TurningPointRight = 0.5 + 1/sqrt(12);
public:
typedef Number value_type;
Number eps, delta;
Number D;
Number dt;
Number uAst;
Number u_regularization;
Number mu;
Number rho_0;
Number rho_1;
Number rho_2;
Number M;
Number M_reduced;
Number fscaling;
Number SurfaceTensionScaling;
//Number rho (const Number& x) {return rho_0;}
Number rho (const Number& x){return x*rho_2 + (1-x)*rho_1;}
Number Sigma1,Sigma2,Sigma3,SigmaT;
// Constructor takes parameter
Problem (const Number& dt_, const Dune::ParameterTree& ptree)
: t(0.0)
{
eps = ptree.get("problem.eps",(Number)0.1);
delta = ptree.get("problem.delta",(Number)0.1);
M = ptree.get("problem.Mobility",(Number)0.1);
M_reduced = ptree.get("problem.MobilityReduced",(Number)0.1);
u_regularization = ptree.get("problem.URegularization",(Number)0.01);
D = ptree.get("physics.D",(Number)0.1);
uAst = ptree.get("physics.uAst",(Number)0.1);
rho_0 = ptree.get("physics.rho",(Number)1);
mu = ptree.get("physics.mu",(Number)0.001);
fscaling = ptree.get("phsyics.fscaling",(Number)1);
SurfaceTensionScaling = ptree.get("physics.SurfaceTensionScaling",(Number)1);
rho_1 = ptree.get("3p.rho1",(Number)1);
rho_2 = ptree.get("3p.rho2",(Number)1);
Sigma1 = ptree.get("3p.Sigma1",(Number)1);
Sigma2 = ptree.get("3p.Sigma2",(Number)1);
Sigma3 = ptree.get("3p.Sigma3",(Number)1);
SigmaT = 1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3);
dt=dt_;
}
Number Mob (const Number& u) const
{
return M; //M_reduced + (M-M_reduced)*q(u);
}
//Double Well
Number W (const Number& u) const
{
return u*u*(1-u)*(1-u);
}
// Double Well Derivative(/Secant)
Number dW (const Number& x, const Number& y) const
{
return (x+y)*(1 +x*x +y*y) -2*(x*x +x*y +y*y);
}
// Double Well Derivative
Number dW (const Number& x) const
{
return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
}
// Double Well Derivative of Convex Part
Number dWConvex (const Number& x) const
{
if(x<TurningPointLeft)
return dW(x)-dW(TurningPointLeft);
if(x>TurningPointRight)
return dW(x)-dW(TurningPointRight);
return 0;
}
// Double Well Derivative of Convcarve Part
Number dWConcarve (const Number& x) const
{
if(x<TurningPointLeft)
return dW(TurningPointLeft);
if(x>TurningPointRight)
return dW(TurningPointRight);
return dW(x);
}
Number dWdpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma1-SigmaT)*dW(pf1)-SigmaT*dW(pf2)-SigmaT*dW(pf3);
}
Number dWdpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma2-SigmaT)*dW(pf2)-SigmaT*dW(pf1)-SigmaT*dW(pf3);
}
Number dWdpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
{
return (Sigma3-SigmaT)*dW(pf3)-SigmaT*dW(pf1)-SigmaT*dW(pf2);
}
Number f (const Number& u) const
{
if(u<0)
return fscaling*(-1.0+1000*u);
else
return fscaling*(u*u-1);
}
Number q (const Number& pf) const
{
if (pf<0 || pf>1)
return 0;
return pf*(1-pf);
}
Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
{
if (pf1<0 || pf3<0)
return 0;
return pf1*pf3;
}
Number Diffusion(const Number& pf_c)
{
return D*fmax(1e-3,(pf_c-0.05)/0.95);
}
Number VelDissipation(const Number& pf_f)
{
return 10*(1-pf_f)*(1-pf_f);
}
//! Set time in instationary case
void setTime (Number t_)
{
t = t_;
}
};
#pragma once
// C includes
#include<sys/stat.h>
// C++ includes
#include<math.h>
#include<iostream>
// dune-common includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
// dune-geometry includes
#include<dune/geometry/referenceelements.hh>
#include<dune/geometry/quadraturerules.hh>
// dune-grid includes
#include<dune/grid/onedgrid.hh>
#include<dune/grid/yaspgrid.hh>
#include<dune/grid/utility/structuredgridfactory.hh>
#include<dune/grid/io/file/vtk.hh>
#include<dune/grid/io/file/gmshreader.hh>
#if HAVE_UG
#include<dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#include<dune/alugrid/dgf.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif
// dune-istl included by pdelab
// dune-pdelab includes
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include<dune/pdelab/finiteelementmap/qkfem.hh>
#include<dune/pdelab/finiteelementmap/p0fem.hh>
#include<dune/pdelab/constraints/common/constraints.hh>
#include<dune/pdelab/constraints/common/constraintsparameters.hh>
#include<dune/pdelab/constraints/conforming.hh>
#include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/interpolate.hh>
#include<dune/pdelab/gridfunctionspace/vtk.hh>
#include<dune/pdelab/gridoperator/gridoperator.hh>
#include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/variablefactories.hh>
#include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/stationary/linearproblem.hh>
#include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/newton/newton.hh>
#include<dune/pdelab/adaptivity/adaptivity.hh>
//#include <dune/xt/grid/grids.hh>
//#include <dune/xt/grid/view/periodic.hh>
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment