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

More cleanup

parent 88a1be34
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune/phasefield/base_parameters.hh
+ 0
− 26
View file @ f7b9e9d4
......@@ -369,29 +369,3 @@ public:
{
}
};
template <typename Number, typename DW, typename Mobility, typename Reaction, typename VelDissipation>
class Params_summerschool : public Params_fs_r<Number, DW, Mobility, Reaction, VelDissipation>
{
public:
Number CellHeight;
Params_summerschool(const Dune::ParameterTree& ptree) :
Params_fs_r<Number, DW, Mobility, Reaction, VelDissipation>(ptree),
CellHeight(ptree.get("domain.CellHeight",0.1))
{
}
};
template <typename Number, typename TW, typename Mobility, typename Reaction, typename VelDissipation>
class Params_summerschool_3p : public Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>
{
public:
Number CellHeight;
Params_summerschool_3p(const Dune::ParameterTree& ptree) :
Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
CellHeight(ptree.get("domain.CellHeight",0.1))
{
}
};
dune/phasefield/boundary_fn.hh deleted 100644 → 0
+ 0
− 79
View file @ 88a1be34
#pragma once
#include <cmath>
#include <dune/pdelab/common/function.hh>
#include <dune/pdelab/constraints/common/constraintsparameters.hh>
#define GRID_EPS 1e-7
#define GRID_WIDTH 2 // 1.5 for summer school //2 for contactpointtest //6 for full summer school
#define GRID_BOTTOM 0 //0 For summer school //0.35 for contactpointtest //0.35 for hainesjump
#define GRID_TOP 1 //0.125 for summer school //1 for contactpointtest //0.65 for hainesjump
#define GRID_TEST 0.05
class InflowOutflow_v : 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]> GRID_WIDTH-GRID_EPS || xg[0] < GRID_EPS) && xg[1]>GRID_BOTTOM+GRID_EPS+GRID_TEST && xg[1]<GRID_TOP-GRID_EPS - GRID_TEST)
{
return false;
}
return true;
}
};
class InflowOutflow_p : 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]> GRID_WIDTH-GRID_EPS || xg[0] < GRID_EPS) && xg[1]>GRID_BOTTOM+GRID_EPS && xg[1]<GRID_TOP-GRID_EPS)
{
return true;
}
return false;
}
};
class InflowOutflow_u : 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]< GRID_EPS)
{
return true;
}
return false;
}
};
class InflowOutflow_phi : 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]< GRID_EPS || xg[0] >GRID_WIDTH-GRID_EPS) && xg[1]>GRID_BOTTOM+GRID_EPS && xg[1]<GRID_TOP-GRID_EPS)
{
return true;
}
return false;
}
};
dune/phasefield/initial_fn.hh
+ 7
− 446
View file @ f7b9e9d4
......@@ -6,6 +6,9 @@
#include <dune/phasefield/base_parameters.hh>
//______________________________________________________________________________________________________________________
template<typename GV, typename RF, std::size_t dim_range, typename Params>
class ZeroInitialVector :
public Dune::PDELab::AnalyticGridFunctionBase<
......@@ -284,316 +287,8 @@ public:
}
};
//______________________________________________________________________________________________________________________
template<typename GV, typename RF, typename Params>
class PhiInitial_2Precip :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_2Precip<GV,RF,Params> >
{
RF eps;
RF x1center;
RF y1center;
RF radius1;
RF x2center;
RF y2center;
RF radius2;
RF x3center;
RF y3center;
RF radius3;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PhiInitial_2Precip(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_2Precip<GV,RF,Params> > (gv)
{
eps = params.eps;
radius1 = params.initial.ptree.get("Radius1",0.2);
x1center = params.initial.ptree.get("x1Center",0.5);
y1center = params.initial.ptree.get("y1Center",0.5);
radius2 = params.initial.ptree.get("Radius2",0.2);
x2center = params.initial.ptree.get("x2Center",0.5);
y2center = params.initial.ptree.get("y2Center",0.5);
radius3 = params.initial.ptree.get("Radius3",0.2);
x3center = params.initial.ptree.get("x3Center",0.5);
y3center = params.initial.ptree.get("y3Center",0.5);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
RF dist1 = std::sqrt( (x[0]-x1center)*(x[0]-x1center) + (x[1]-y1center)*(x[1]-y1center) );
RF y1 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist1-radius1)));
RF dist2 = std::sqrt( (x[0]-x2center)*(x[0]-x2center) + (x[1]-y2center)*(x[1]-y2center) );
RF y2 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist2-radius2)));
RF dist3 = std::sqrt( (x[0]-x3center)*(x[0]-x3center) + (x[1]-y3center)*(x[1]-y3center) );
RF y3 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist3-radius3)));
y[0] = y1*y2*y3;
}
};
template<typename GV, typename RF, typename Params>
class PhiInitial_SummerSchoolFull :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_SummerSchoolFull<GV,RF,Params> >
{
public:
RF eps;
RF radius;
std::vector<RF> xValues;
std::vector<RF> yValues;
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PhiInitial_SummerSchoolFull(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_SummerSchoolFull<GV,RF,Params> > (gv)
{
eps = params.eps;
radius = params.initial.ptree.get("Radius",0.2);
xValues.push_back(0.560); yValues.push_back(0.070);
xValues.push_back(0.831); yValues.push_back(0.029);
xValues.push_back(1.353); yValues.push_back(0.008);
xValues.push_back(1.298); yValues.push_back(0.105);
xValues.push_back(1.387); yValues.push_back(0.108);
xValues.push_back(1.959); yValues.push_back(0.014);
xValues.push_back(2.397); yValues.push_back(-0.001);
xValues.push_back(3.014); yValues.push_back(0.074);
xValues.push_back(3.703); yValues.push_back(0.244);
xValues.push_back(3.808); yValues.push_back(-0.015);
xValues.push_back(4.190); yValues.push_back(0.071);
xValues.push_back(4.816); yValues.push_back(0.026);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
RF dist;
RF temp;
RF yc = 1;
for(int i=0; i<xValues.size(); i++)
{
dist = std::sqrt( (x[0]-xValues[i])*(x[0]-xValues[i]) + (x[1]-yValues[i])*(x[1]-yValues[i]) );
temp = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-radius)));
yc = yc*temp;
}
y[0] = yc;
}
};
template<typename GV, typename RF, typename Params, const int phase = -1>
class PhiInitial_SummerSchoolPaper :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_SummerSchoolPaper<GV,RF,Params,phase> >
{
public:
RF eps;
RF radius;
std::vector<RF> xValues;
std::vector<RF> yValues;
private:
inline static RF PxToCoord_x(const RF px_x)
{
return 0.5+(px_x - 150)/((RF)299);
//return 0.5+(px_x - 150)/((RF)299)-3.5; //-> Christians Vortrag
}
inline static RF PxToCoord_y(const RF px_y)
{
return (102 - px_y)/((RF)299);
}
inline static RF PxToCoord_x2(const RF px_x)
{
return (px_x - 163)/((RF)314);
}
inline static RF PxToCoord_y2(const RF px_y)
{
return (107 - px_y)/((RF)314);
}
inline static RF PxToCoord_x3(const RF px_x)
{
return (px_x - 391)/((RF)224/1.5);
}
inline static RF PxToCoord_y3(const RF px_y)
{
return (749 - px_y)/((RF)224/1.5);
}
inline static RF PxToCoord_x4(const RF px_x)
{
return (px_x - 1081)/((RF)189/0.75);
}
inline static RF PxToCoord_y4(const RF px_y)
{
return (2037 - px_y)/((RF)189/0.75);
}
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PhiInitial_SummerSchoolPaper(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_SummerSchoolPaper<GV,RF,Params,phase> > (gv)
{
eps = params.eps;
radius = params.initial.ptree.get("Radius",0.2);
//Experiment 1
//xValues.push_back(PxToCoord_x(114)); yValues.push_back(PxToCoord_y(81)); //Extra Nucleus, to check the influence on Nucleus 1
xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75)); //Nucleus 1
xValues.push_back(PxToCoord_x(317)); yValues.push_back(PxToCoord_y(85));
xValues.push_back(PxToCoord_x(364)); yValues.push_back(PxToCoord_y(110));
xValues.push_back(PxToCoord_x(556)); yValues.push_back(PxToCoord_y(65));
xValues.push_back(PxToCoord_x(759)); yValues.push_back(PxToCoord_y(85));
xValues.push_back(PxToCoord_x(832)); yValues.push_back(PxToCoord_y(159));
xValues.push_back(PxToCoord_x(938)); yValues.push_back(PxToCoord_y(101));
xValues.push_back(PxToCoord_x(1100)); yValues.push_back(PxToCoord_y(62));
// xValues.push_back(PxToCoord_x(1240)); yValues.push_back(PxToCoord_y(70)); //Nuclei für Christians Vortrag
xValues.push_back(PxToCoord_x(1229)); yValues.push_back(PxToCoord_y(20));
xValues.push_back(PxToCoord_x(1258)); yValues.push_back(PxToCoord_y(29));
xValues.push_back(PxToCoord_x(1340)); yValues.push_back(PxToCoord_y(57));
xValues.push_back(PxToCoord_x(1346)); yValues.push_back(PxToCoord_y(77));
xValues.push_back(PxToCoord_x(1354)); yValues.push_back(PxToCoord_y(103));
xValues.push_back(PxToCoord_x(1731)); yValues.push_back(PxToCoord_y(115));
xValues.push_back(PxToCoord_x(1915)); yValues.push_back(PxToCoord_y(77));
// Addional Nuclei
//xValues.push_back(PxToCoord_x(668)); yValues.push_back(PxToCoord_y(63));
//xValues.push_back(PxToCoord_x(775)); yValues.push_back(PxToCoord_y(143));
//xValues.push_back(PxToCoord_x(874)); yValues.push_back(PxToCoord_y(139));
//xValues.push_back(PxToCoord_x(815)); yValues.push_back(PxToCoord_y(36));
//xValues.push_back(PxToCoord_x(848)); yValues.push_back(PxToCoord_y(53));
//xValues.push_back(PxToCoord_x(1432)); yValues.push_back(PxToCoord_y(101));
/*
//Experiment 7
xValues.push_back(PxToCoord_x2(423)); yValues.push_back(PxToCoord_y2(129));
xValues.push_back(PxToCoord_x2(1636)); yValues.push_back(PxToCoord_y2(98));
*/
/*
//Inflow
xValues.push_back(PxToCoord_x3(402)); yValues.push_back(PxToCoord_y3(642));
xValues.push_back(PxToCoord_x3(406)); yValues.push_back(PxToCoord_y3(689));
xValues.push_back(PxToCoord_x3(416)); yValues.push_back(PxToCoord_y3(696));
xValues.push_back(PxToCoord_x3(425)); yValues.push_back(PxToCoord_y3(709));
xValues.push_back(PxToCoord_x3(441)); yValues.push_back(PxToCoord_y3(709));
xValues.push_back(PxToCoord_x3(468)); yValues.push_back(PxToCoord_y3(719));
xValues.push_back(PxToCoord_x3(435)); yValues.push_back(PxToCoord_y3(681));
xValues.push_back(PxToCoord_x3(460)); yValues.push_back(PxToCoord_y3(673));
xValues.push_back(PxToCoord_x3(476)); yValues.push_back(PxToCoord_y3(672));
xValues.push_back(PxToCoord_x3(497)); yValues.push_back(PxToCoord_y3(693));
xValues.push_back(PxToCoord_x3(548)); yValues.push_back(PxToCoord_y3(687));
xValues.push_back(PxToCoord_x3(594)); yValues.push_back(PxToCoord_y3(680));
*/
/*
//Inflow S1
xValues.push_back(PxToCoord_x4(1090)); yValues.push_back(PxToCoord_y4(1858));
xValues.push_back(PxToCoord_x4(1102)); yValues.push_back(PxToCoord_y4(1885));
xValues.push_back(PxToCoord_x4(1126)); yValues.push_back(PxToCoord_y4(1892));
xValues.push_back(PxToCoord_x4(1134)); yValues.push_back(PxToCoord_y4(1905));
xValues.push_back(PxToCoord_x4(1186)); yValues.push_back(PxToCoord_y4(1901));
xValues.push_back(PxToCoord_x4(1248)); yValues.push_back(PxToCoord_y4(1899));
xValues.push_back(PxToCoord_x4(1296)); yValues.push_back(PxToCoord_y4(1908));
xValues.push_back(PxToCoord_x4(1093)); yValues.push_back(PxToCoord_y4(1951));
xValues.push_back(PxToCoord_x4(1122)); yValues.push_back(PxToCoord_y4(1936));
xValues.push_back(PxToCoord_x4(1143)); yValues.push_back(PxToCoord_y4(1942));
xValues.push_back(PxToCoord_x4(1095)); yValues.push_back(PxToCoord_y4(1975));
xValues.push_back(PxToCoord_x4(1141)); yValues.push_back(PxToCoord_y4(1970));
xValues.push_back(PxToCoord_x4(1155)); yValues.push_back(PxToCoord_y4(1985));
xValues.push_back(PxToCoord_x4(1171)); yValues.push_back(PxToCoord_y4(1967));
xValues.push_back(PxToCoord_x4(1230)); yValues.push_back(PxToCoord_y4(1949));
xValues.push_back(PxToCoord_x4(1249)); yValues.push_back(PxToCoord_y4(1957));
xValues.push_back(PxToCoord_x4(1273)); yValues.push_back(PxToCoord_y4(1963));
xValues.push_back(PxToCoord_x4(1122)); yValues.push_back(PxToCoord_y4(2018));
xValues.push_back(PxToCoord_x4(1134)); yValues.push_back(PxToCoord_y4(2010));
xValues.push_back(PxToCoord_x4(1184)); yValues.push_back(PxToCoord_y4(2021));
xValues.push_back(PxToCoord_x4(1199)); yValues.push_back(PxToCoord_y4(1998));
xValues.push_back(PxToCoord_x4(1283)); yValues.push_back(PxToCoord_y4(2029));
*/
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
RF dist;
RF temp;
RF yc = 1;
for(int i=0; i<xValues.size(); i++)
{
dist = std::sqrt( (x[0]-xValues[i])*(x[0]-xValues[i]) + (x[1]-yValues[i])*(x[1]-yValues[i]) );
temp = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-radius)));
yc = yc*temp;
}
if(phase >= 0)
{
RF yf = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-0.3)));
if (phase == 0)
{
yc = yc*yf;
}
else
{
yc = yc*(1-yf);
}
}
y[0] = yc;
}
};
template<typename GV, typename RF, typename Params, const int Phase>
class PhiInitial_HainesJump :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_HainesJump<GV,RF,Params,Phase> >
{
RF eps;
RF period;
RF amplitude;
RF height;
RF offset;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PhiInitial_HainesJump(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_HainesJump<GV,RF,Params,Phase> > (gv)
{
eps = params.eps;
period = params.initial.ptree.get("Period",1);
amplitude = params.initial.ptree.get("Amplitude",0.3);
height = params.initial.ptree.get("Height",1);
offset = params.initial.ptree.get("Offset",0.1);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
RF temp = (-std::cos(x[0]*2*3.1415*period) + 1)*amplitude/2 + offset;
RF fluid = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-temp)));
fluid = fluid*0.5*(1+tanh(1/eps/sqrt(2)*(height-x[1]-temp)));
RF d1 = x[0]-0.25;
if(Phase == 0)
d1 = -d1;
y=0.5*(1+tanh(1/eps/sqrt(2)*d1))*fluid;
}
};
//______________________________________________________________________________________________________________________
template<typename GV, typename RF, typename Params>
class PhiInitial_SpinodalDecomposition :
......@@ -622,6 +317,9 @@ public:
}
};
//______________________________________________________________________________________________________________________
template<typename GV, typename RF, typename Params, const int Phase>
class PhiInitial_3pBubble :
public Dune::PDELab::AnalyticGridFunctionBase<
......@@ -655,140 +353,3 @@ public:
y=fluid*rf;
}
};
//______________________________________________________________________________________________________________________
#include <dune/phasefield/boundary_fn.hh>
template<typename GV, typename RF, typename Params, int dim>
class VelocityInitial_Inflow :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
VelocityInitial_Inflow<GV,RF,Params,dim> >
{
const Params_Time& time;
RF meanflow;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
VelocityInitial_Inflow(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,VelocityInitial_Inflow<GV,RF,Params,dim> > (gv)
, time(params.time)
{
meanflow = params.initial.ptree.get("MeanInflow",(RF)0.1);
}
template <typename T>
void setTime(T t){
time = t;
}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
if(x[0]<GRID_EPS) //Inflow boundary
{
RF factor = meanflow * 6/( (GRID_TOP-GRID_BOTTOM)*(GRID_TOP-GRID_BOTTOM) );
y[0] = std::max((RF)0, factor * (x[1]-GRID_BOTTOM)*(GRID_TOP-x[1]));
}
}
};
template<typename GV, typename RF, typename Params, int dim>
class VelocityInitial_Channelflow :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
VelocityInitial_Channelflow<GV,RF,Params,dim> >
{
const Params_Time& time;
RF v;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
typedef typename Traits::DomainType DomainType;
typedef typename Traits::RangeType RangeType;
VelocityInitial_Channelflow(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,VelocityInitial_Channelflow<GV,RF,Params,dim> > (gv)
, time(params.time)
{
v = params.initial.ptree.get("FlowV",(RF)0.1);
}
template <typename T>
void setTime(T t){
time = t;
}
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
y = 0.0;
if(x[1]>1-GRID_EPS) //Flow boundary
{
y[0] = v;
}
}
};
template<typename GV, typename RF, typename Params>
class UInitial_Inflow :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, UInitial_Inflow<GV,RF,Params> >
{
RF uConst;
RF uInflow;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
UInitial_Inflow(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,UInitial_Inflow<GV,RF,Params> > (gv)
{
uConst = params.initial.ptree.get("uConst",(RF)0.5);
uInflow = params.initial.ptree.get("uInflow",(RF)0.5);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
y = uConst;
if(x[0]<1e-7)
{
y=uInflow;
}
}
};
template<typename GV, typename RF, typename Params>
class PInitial_Inflow :
public Dune::PDELab::AnalyticGridFunctionBase<
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PInitial_Inflow<GV,RF,Params> >
{
RF pConst;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PInitial_Inflow(const GV & gv, const Params & params) :
Dune::PDELab::AnalyticGridFunctionBase<Traits,PInitial_Inflow<GV,RF,Params> > (gv)
{
pConst = params.initial.ptree.get("pConst",(RF)0.01);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
y = 0;
if(x[0]<1e-7)
{
y=pConst;
}
}
};
dune/phasefield/localoperator/fem_fs_chns_u_only.hh deleted 100644 → 0
+ 0
− 110
View file @ 88a1be34
#pragma once
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/finiteelement/localbasiscache.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/gridfunctionspace/gridfunctionspace.hh>
#include<dune/phasefield/localoperator/fem_base_chns.hh>
//Parameter class (for functions), FEM space and RangeField
template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
class FEM_FS_CHNS_U_Only :
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_U_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_U_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
{
private:
typedef typename CHContainer::field_type RF;
typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
CHLocalBasisCache<CHGFS> chCache;
VLocalBasisCache<NSGFS> vCache;
typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
mutable CHLView chOldLocal;
typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
mutable NSLView nsLocal;
Param& param; // parameter functions
public:
enum {doPatternVolume=true};
enum {doAlphaVolume=true};
//constructor
FEM_FS_CHNS_U_Only(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
chOldLocal(chGfs_,chVOld_),
nsLocal(nsGfs_,nsV),
param(param_)
{
}
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
const int order = 5; //TODO: Choose more educated
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
//chOldLocal.DoEvaluation(eg);
nsLocal.DoSubEvaluation(eg,child(nsLocal.lfs,Dune::Indices::_0));
for(const auto& ip : rule)
{
const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
CHVariables_U_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(),chOldLocal.lv, S);
NSVariables<RF,dim,decltype(nsLocal.lfs)> ns(vCache, ip.position(), nsLocal.lfs, nsLocal.lv, S);
RF pf_c = ch.pf;
RF pf_f = ch.pf + 2*param.delta*(1-ch.pf);
Dune::FieldVector<RF,dim> gradpf_f(0.0);
gradpf_f.axpy((1-2*param.delta),ch.gradpf);
RF pf_c_Old = ch.pfOld;
RF pf_c_reg = pf_c+param.delta;
RF pf_c_reg_Old = pf_c_Old + param.delta;
for (size_t i=0; i<ch.pfspace.size(); i++)
{
RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
RF R1 = -Rspeed/param.eps*
(param.reaction.eval(ch.u)+param.phys.curvatureAlpha*ch.xi)*ch.basis.phi[i]; //TODO Scale with sigma13
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
r.accumulate(ch.pfspace,i,factor*(
(ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
));
r.accumulate(ch.xispace,i,factor*(
(ch.xi-0)*ch.basis.phi[i]
));
// integrate (pf-pfold)*(u-uast)*phi_i/delta_t + pfold*(u-uold)*phi_i/delta_t + D*pf*gradu*grad phi_i = 0
r.accumulate(ch.uspace,i,factor*(
//(u-param.uAst)*phi[i]/param.dt //For Testing
R1 *(ch.u-param.phys.uAst)
+ pf_c_reg_Old*(ch.u-ch.uOld)*ch.basis.phi[i]/param.time.dt
+ param.phys.D*pf_c_reg_Old*(ch.gradu*ch.basis.gradphi[i][0])
+ ( pf_f*(ns.v*ch.gradu) )*ch.basis.phi[i]
- pf_c_reg*param.reaction.eval_bulkreaction()*ch.basis.phi[i]//Bulk reaction
));
}
}
}
};
dune/phasefield/localoperator/fem_fs_r_only.hh deleted 100644 → 0
+ 0
− 83
View file @ 88a1be34
#pragma once
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/finiteelement/localbasiscache.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/gridfunctionspace/gridfunctionspace.hh>
#include<dune/phasefield/localoperator/fem_base_chns.hh>
//Parameter class (for functions), FEM space and RangeField
template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
class FEM_FS_R_Only :
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
public Dune::PDELab::NumericalJacobianVolume<FEM_FS_R_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_R_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
{
private:
typedef typename CHContainer::field_type RF;
typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
CHLocalBasisCache<CHGFS> chCache;
typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
mutable CHLView chOldLocal;
Param& param; // parameter functions
public:
enum {doPatternVolume=true};
enum {doAlphaVolume=true};
//constructor
FEM_FS_R_Only(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
chOldLocal(chGfs_,chVOld_),
param(param_)
{
}
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
const int order = 5; //TODO: Choose more educated
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
for(const auto& ip : rule)
{
const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
CHVariables_U_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(),chOldLocal.lv, S);
for (size_t i=0; i<ch.pfspace.size(); i++)
{
RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
RF R1 = -Rspeed/param.eps*
(param.reaction.eval(ch.uOld))*ch.basis.phi[i];
r.accumulate(ch.pfspace,i,factor*(
(ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
-R1
));
r.accumulate(ch.xispace,i,factor*(
(ch.xi-0)*ch.basis.phi[i]
));
r.accumulate(ch.uspace,i,factor*(
(ch.u-ch.uOld)*ch.basis.phi[i]
));
}
}
}
};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment