diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh index 6d8b0a0bc1708be66da1051234eadf2361b6e2cf..a63a2031dba920d32030613ebfe810a640a637b9 100644 --- a/dune/phasefield/base_parameters.hh +++ b/dune/phasefield/base_parameters.hh @@ -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)) - { - } -}; diff --git a/dune/phasefield/boundary_fn.hh b/dune/phasefield/boundary_fn.hh deleted file mode 100644 index b6592c77dab41aeea4ca460fc90f338073d6d544..0000000000000000000000000000000000000000 --- a/dune/phasefield/boundary_fn.hh +++ /dev/null @@ -1,79 +0,0 @@ - -#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; - } -}; diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh index 8b98df41364bf4f2a4fb045d6f3e651003977690..04488b71270a938435544d23e8576ae2afb2563d 100644 --- a/dune/phasefield/initial_fn.hh +++ b/dune/phasefield/initial_fn.hh @@ -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; - } - } -}; diff --git a/dune/phasefield/localoperator/fem_fs_chns_u_only.hh b/dune/phasefield/localoperator/fem_fs_chns_u_only.hh deleted file mode 100644 index c1fca5a1f23068004682037b02ad9cef123ed6ba..0000000000000000000000000000000000000000 --- a/dune/phasefield/localoperator/fem_fs_chns_u_only.hh +++ /dev/null @@ -1,110 +0,0 @@ -#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 - )); - } - } - } -}; diff --git a/dune/phasefield/localoperator/fem_fs_r_only.hh b/dune/phasefield/localoperator/fem_fs_r_only.hh deleted file mode 100644 index 902946f1460be96f6eaaff87f46a9169bffeca37..0000000000000000000000000000000000000000 --- a/dune/phasefield/localoperator/fem_fs_r_only.hh +++ /dev/null @@ -1,83 +0,0 @@ -#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] - )); - } - } - } -};