diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 2089d4efeac4d16e24d7dea721c48d6aab4eb3fa..e35943de17b28b1dffe3f0cf727f0c5dfaf20f59 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -1,202 +1,13 @@
#pragma once
//#include "../myincludes.hh"
+#include <dune/phasefield/base_potentials.hh>
#include <cmath>
#include <limits>
#include <dune/common/parametertreeparser.hh>
#include <iostream>
-template<typename Number>
-class DoubleWell_poly
-{
-public:
- const Number TurningPointLeft = 0.5 - 1/sqrt(12);
- const Number TurningPointRight = 0.5 + 1/sqrt(12);
-
- DoubleWell_poly(const Dune::ParameterTree & ptree)
- {
- }
-
- //Double Well
- Number eval(const Number& x) const
- {
- return x*x*(1-x)*(1-x);
- }
-
- // Double Well Derivative
- Number eval_d(const Number& x) const
- {
- return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
- }
-
- // Double Well Derivative(/Secant)
- Number eval_d (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 of Convex Part
- Number eval_dConvex (const Number& x) const
- {
- if(x<TurningPointLeft)
- return eval_d(x)- eval_d(TurningPointLeft);
- if(x>TurningPointRight)
- return eval_d(x)-eval_d(TurningPointRight);
- return 0;
- }
-
- // Double Well Derivative of Convcarve Part
- Number eval_dConcarve (const Number& x) const
- {
- if(x<TurningPointLeft)
- return eval_d(TurningPointLeft);
- if(x>TurningPointRight)
- return eval_d(TurningPointRight);
- return eval_d(x);
- }
-
- // The function satisfying q=sqrt(2 W), inbetween 0 and 1
- Number eval_q(const Number& x) const
- {
- if(x<0 || x>1)
- return 0;
- else
- return sqrt(2)*x*(1-x);
- }
-};
-
-
-template<typename Number>
-class DoubleWell_limited
-{
-private:
- DoubleWell_poly<Number> wp;
-
- Number limiter(const Number& x) const
- {
- if(x<=-delta)
- {
- return std::numeric_limits<Number>::max();
- }
- if(x>=0)
- {
- return 0;
- }
- return x*x/(delta+x);
- }
-
- Number limiter_d(const Number& x) const
- {
- if(x<=-delta)
- {
- return std::numeric_limits<Number>::max();
- }
- if(x>=0)
- {
- return 0;
- }
- return x*(x+2*delta)/((x+delta)*(x+delta));
- }
-
-public:
- const Number delta;
-
- DoubleWell_limited(const Dune::ParameterTree & ptree) :
- wp(ptree),
- delta(ptree.get("DoubleWellDelta",(Number)0.1))
- {
- }
-
- //Double Well
- Number eval(const Number& x) const
- {
- return wp.eval(x) + limiter(x) + limiter(1-x);
- }
-
- // Double Well Derivative
- Number eval_d(const Number& x) const
- {
- return wp.eval_d(x) + limiter_d(x) + limiter_d(1-x);
- }
-
- // Double Well Derivative(/Secant)
- Number eval_d (const Number& x, const Number& y) const
- {
- std::cout << "DoubleWell_limited secant not properly implemented yet" << std::endl;
- if(x==y)
- {
- return eval_d(x);
- }
- else
- {
- return (eval(x)-eval(y))/(x-y);
- }
- }
-
- // Double Well Derivative of Convex Part
- Number eval_dConvex (const Number& x) const
- {
- return wp.eval_dConvex(x) + limiter_d(x) + limiter_d(1-x);
- }
-
- // Double Well Derivative of Convcarve Part
- Number eval_dConcarve (const Number& x) const
- {
- return wp.eval_dConcarve(x);
- }
-
- // Just foward q, we do not care about the limiter for this
- Number eval_q (const Number& x) const
- {
- return wp.eval_q(x);
- }
-};
-
-template<typename Number, typename DW>
-class TripleWell
-{
-public:
- const Number Sigma1;
- const Number Sigma2;
- const Number Sigma3;
- const Number SigmaT;
-
- DW dw;
-
- TripleWell(const Dune::ParameterTree & param) :
- Sigma1(param.get("Sigma1",(Number)1)),
- Sigma2(param.get("Sigma2",(Number)1)),
- Sigma3(param.get("Sigma3",(Number)1)),
- SigmaT(1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3)),
- dw(param)
- {
- }
-
- Number eval_dpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
- {
- return (Sigma1-SigmaT)*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf2)-SigmaT*dw.eval_d(pf3);
- }
-
- Number eval_dpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
- {
- return (Sigma2-SigmaT)*dw.eval_d(pf2)-SigmaT*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf3);
- }
-
- Number eval_dpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
- {
- return (Sigma3-SigmaT)*dw.eval_d(pf3)-SigmaT*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf2);
- }
-
- Number eval_q (const Number& pf1, const Number& pf2, const Number& pf3) const
- {
- if(abs(pf2-1)<1e-14)
- return 0;
- return dw.eval_q(pf1/(1-pf2))*(1-pf2)*(1-pf2); //A function being q on the 13 interface and 0 on the others
- }
-};
-
-
//______________________________________________________________________________________________________
template<typename Number>
diff --git a/dune/phasefield/base_potentials.hh b/dune/phasefield/base_potentials.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f170addc1f0f99c79b7846aa934c1595cd55af7b
--- /dev/null
+++ b/dune/phasefield/base_potentials.hh
@@ -0,0 +1,265 @@
+
+#pragma once
+//#include "../myincludes.hh"
+#include <dune/phasefield/base_potentials.hh>
+#include <cmath>
+#include <limits>
+#include <dune/common/parametertreeparser.hh>
+#include <iostream>
+
+
+//______________________________________________________________________________________________________
+
+
+template<typename Number>
+class DoubleWell_poly
+{
+public:
+ const Number TurningPointLeft = 0.5 - 1/sqrt(12);
+ const Number TurningPointRight = 0.5 + 1/sqrt(12);
+
+ DoubleWell_poly(const Dune::ParameterTree & ptree)
+ {
+ }
+
+ //Double Well
+ Number eval(const Number& x) const
+ {
+ return x*x*(1-x)*(1-x);
+ }
+
+ // Double Well Derivative
+ Number eval_d(const Number& x) const
+ {
+ return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
+ }
+
+ // Double Well Derivative(/Secant)
+ Number eval_d (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 of Convex Part
+ Number eval_dConvex (const Number& x) const
+ {
+ if(x<TurningPointLeft)
+ return eval_d(x)- eval_d(TurningPointLeft);
+ if(x>TurningPointRight)
+ return eval_d(x)-eval_d(TurningPointRight);
+ return 0;
+ }
+
+ // Double Well Derivative of Convcarve Part
+ Number eval_dConcarve (const Number& x) const
+ {
+ if(x<TurningPointLeft)
+ return eval_d(TurningPointLeft);
+ if(x>TurningPointRight)
+ return eval_d(TurningPointRight);
+ return eval_d(x);
+ }
+
+ // The function satisfying q=sqrt(2 W), inbetween 0 and 1
+ Number eval_q(const Number& x) const
+ {
+ if(x<0 || x>1)
+ return 0;
+ else
+ return sqrt(2)*x*(1-x);
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template<typename Number>
+class Limiter_Diverging
+{
+public:
+ const Number delta;
+
+
+ Limiter_Diverging(const Dune::ParameterTree & ptree) :
+ delta(ptree.get("DoubleWellDelta",(Number)0.1))
+ {
+ }
+
+ Number eval(const Number& x) const
+ {
+ if(x<=-delta)
+ {
+ return std::numeric_limits<Number>::max();
+ }
+ if(x>=0)
+ {
+ return 0;
+ }
+ return x*x/(delta+x);
+ }
+
+ Number eval_d(const Number& x) const
+ {
+ if(x<=-delta)
+ {
+ return std::numeric_limits<Number>::max();
+ }
+ if(x>=0)
+ {
+ return 0;
+ }
+ return x*(x+2*delta)/((x+delta)*(x+delta));
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template<typename Number>
+class Limiter_FakeDiverging
+{
+public:
+ const Number delta; //Where the potential would be diverging
+ const Number gamma; //Where we go into a linear function
+ Number gammaeval;
+ Number gammaeval_d;
+ Limiter_FakeDiverging(const Dune::ParameterTree & ptree) :
+ delta(ptree.get("DoubleWellDelta",(Number)0.1)),
+ gamma(ptree.get("DoubleWellGamma",(Number)delta/20))
+ {
+ gammaeval=delta*gamma*gamma/(delta-gamma);
+ gammaeval_d=-delta*gamma*(-gamma+2*delta)/((delta-gamma)*(delta-gamma));
+ }
+
+ Number eval(const Number& x) const
+ {
+ if(x<=-gamma)
+ {
+ std::cout << "Used FakeDiverging" <<std::endl;
+ return gammaeval + (x-gamma)*gammaeval_d;
+ }
+ if(x>=0)
+ {
+ return 0;
+ }
+ return delta*x*x/(delta+x);
+ }
+
+ Number eval_d(const Number& x) const
+ {
+ if(x<=-gamma)
+ {
+ //std::cout << "Used FakeDiverging" <<std::endl;
+ return gammaeval_d;
+ }
+ if(x>=0)
+ {
+ return 0;
+ }
+ return delta*x*(x+2*delta)/((x+delta)*(x+delta));
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template<typename Number, typename DoubleWell, typename Limiter>
+class DoubleWell_limited
+{
+public:
+ DoubleWell wp;
+ Limiter limiter;
+
+ DoubleWell_limited(const Dune::ParameterTree & ptree) :
+ wp(ptree),
+ limiter(ptree)
+ {
+ }
+
+ //Double Well
+ Number eval(const Number& x) const
+ {
+ return wp.eval(x) + limiter.eval(x) + limiter.eval(1-x);
+ }
+
+ // Double Well Derivative
+ Number eval_d(const Number& x) const
+ {
+ return wp.eval_d(x) + limiter.eval_d(x) - limiter.eval_d(1-x);
+ }
+
+ // Double Well Derivative(/Secant)
+ Number eval_d (const Number& x, const Number& y) const
+ {
+ std::cout << "DoubleWell_limited secant not properly implemented yet" << std::endl;
+ if(x==y)
+ {
+ return eval_d(x);
+ }
+ else
+ {
+ return (eval(x)-eval(y))/(x-y);
+ }
+ }
+
+ // Double Well Derivative of Convex Part
+ Number eval_dConvex (const Number& x) const
+ {
+ return wp.eval_dConvex(x) + limiter.eval_d(x) - limiter.eval_d(1-x);
+ }
+
+ // Double Well Derivative of Convcarve Part
+ Number eval_dConcarve (const Number& x) const
+ {
+ return wp.eval_dConcarve(x);
+ }
+
+ // Just foward q, we do not care about the limiter for this
+ Number eval_q (const Number& x) const
+ {
+ return wp.eval_q(x);
+ }
+};
+
+
+//______________________________________________________________________________________________________
+
+template<typename Number, typename DW>
+class TripleWell
+{
+public:
+ const Number Sigma1;
+ const Number Sigma2;
+ const Number Sigma3;
+ const Number SigmaT;
+
+ DW dw;
+
+ TripleWell(const Dune::ParameterTree & param) :
+ Sigma1(param.get("Sigma1",(Number)1)),
+ Sigma2(param.get("Sigma2",(Number)1)),
+ Sigma3(param.get("Sigma3",(Number)1)),
+ SigmaT(1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3)),
+ dw(param)
+ {
+ }
+
+ Number eval_dpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma1-SigmaT)*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf2)-SigmaT*dw.eval_d(pf3);
+ }
+
+ Number eval_dpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma2-SigmaT)*dw.eval_d(pf2)-SigmaT*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf3);
+ }
+
+ Number eval_dpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma3-SigmaT)*dw.eval_d(pf3)-SigmaT*dw.eval_d(pf1)-SigmaT*dw.eval_d(pf2);
+ }
+
+ Number eval_q (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ if(abs(pf2-1)<1e-14)
+ return 0;
+ return dw.eval_q(pf1/(1-pf2))*(1-pf2)*(1-pf2); //A function being q on the 13 interface and 0 on the others
+ }
+};
diff --git a/dune/phasefield/boundary_fn.hh b/dune/phasefield/boundary_fn.hh
index 06ea10116fa05cf11172cb05916a963076def580..becbe5700cac6dd7195195038018a4e902d28463 100644
--- a/dune/phasefield/boundary_fn.hh
+++ b/dune/phasefield/boundary_fn.hh
@@ -22,6 +22,8 @@ public:
}
};
+
+
class InflowOutflow_p : public Dune::PDELab::DirichletConstraintsParameters
{
public:
diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index 2405e7dd4c5a1d359e08110247296a917bdc1f4b..38d762b67d704b9a62a8d0f9b0e037f78c2473f0 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -69,7 +69,7 @@ public:
{
if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
ns_newton.apply();
-
+
if (verb > 1) std::cout << "= Begin newtonapply CH:" << std::endl;
ch_newton.apply();
@@ -114,6 +114,8 @@ public:
}
};
+//_________________________________________________________________________________________________________________________
+
template<typename GFS, typename NS_BDRY_FN, typename CH_BDRY_FN, typename NS_V, typename CH_V>
void interpolateToBdry(const GFS& gfs, const NS_BDRY_FN& nsBdryFn, const CH_BDRY_FN& chBdryFn, NS_V& nsV, CH_V& chV)
{
diff --git a/dune/phasefield/fs_chns_r.hh b/dune/phasefield/fs_chns_r.hh
index f88ad39718b4954eb33de720c32600aad4419899..e2e04a39f17544c9bfed72addc0e9dcee5668d66 100644
--- a/dune/phasefield/fs_chns_r.hh
+++ b/dune/phasefield/fs_chns_r.hh
@@ -196,3 +196,36 @@ public:
{
}
};
+
+template<typename BdryV1, typename BdryV2, typename BdryP, typename BdryPhi, typename BdryMu, typename BdryU>
+class Bdry_fs_chns_r_compositeV_2d
+{
+public:
+ const int dim=2;
+
+ BdryV1 bdryV1;
+ BdryV2 bdryV2;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ BdryU bdryU;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu, BdryU> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ Bdry_fs_chns_r_compositeV_2d() :
+ bdryV1(),
+ bdryV2(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryU(),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryU)
+ {
+ }
+};
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index 265d77f3c86230ed040fe53ab9490afd4a0b9105..a9522f69f7b841857d95e4fee56890023c3ea0e7 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -146,6 +146,36 @@ public:
//______________________________________________________________________________________________________________________
+template<typename GV, typename RF, typename Params>
+class PhiInitial_Rectangle :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Rectangle<GV,RF,Params> >
+{
+ RF eps;
+ RF radius;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_Rectangle(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Rectangle<GV,RF,Params> > (gv)
+ {
+ eps = params.eps;
+ radius = params.initial.ptree.get("Radius",0.2);
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ RF dist1 = std::abs(x[0]-0.5);
+ RF dist2 = std::abs(x[1]-0.5);
+ RF m1 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist1-radius)));
+ RF m2 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist2-radius)));
+ y[0] = 1-(1-m1)*(1-m2);
+ }
+};
+
+//______________________________________________________________________________________________________________________
+
template<typename GV, typename RF, typename Params, int dim>
class VelocityInitial_Inflow :
public Dune::PDELab::AnalyticGridFunctionBase<
@@ -179,7 +209,7 @@ public:
y = 0.0;
if(x[0]<1e-7) //Inflow boundary
{
- y[0] = std::max((RF)0,3/2 * meanflow * (x[1]+1)*(1-x[1]));
+ y[0] = std::max((RF)0,3/2 * meanflow * x[1]*(1-x[1]));
}
}
};
diff --git a/dune/phasefield/localoperator/fem_base_chns.hh b/dune/phasefield/localoperator/fem_base_chns.hh
index a6ef66e7630a2bf595d56fdf541974ad5aaa0ea1..3826fbcd173f7d495b98604bb96b34da15c0766c 100644
--- a/dune/phasefield/localoperator/fem_base_chns.hh
+++ b/dune/phasefield/localoperator/fem_base_chns.hh
@@ -5,10 +5,10 @@
#include<dune/pdelab/localoperator/variablefactories.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
-
template<typename GFS, typename Container, typename Entity, typename RangeFieldType>
class GFS_LocalViewer
{
+public:
typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
mutable LFS lfs;
typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
@@ -18,21 +18,21 @@ class GFS_LocalViewer
mutable std::vector<RangeFieldType> lv;
- Entity& entity;
+ Entity entity;
GFS_LocalViewer(GFS& gfs, Container& v) :
lfs(gfs), lfsCache(lfs), lView(v), lv(gfs.maxLocalSize())
{
}
- template<typename Femspace>
- bool DoSubEvaluation(const Entity& newEntity, const Femspace& femspace)
+ template<typename EG, typename Femspace>
+ bool DoSubEvaluation(const EG& eg, const Femspace& femspace)
{
- if(entity == newEntity)
+ if(entity == eg.entity())
{
return false;
}
- entity = newEntity;
+ entity = eg.entity();
lfs.bind(entity);
lfsCache.update();
@@ -44,13 +44,14 @@ class GFS_LocalViewer
return true;
}
- bool DoEvaluation(const Entity& newEntity)
+ template<typename EG>
+ bool DoEvaluation(const EG& eg)
{
- if(entity == newEntity)
+ if(entity == eg.entity())
{
return false;
}
- entity = newEntity;
+ entity = eg.entity();
lfs.bind(entity);
lfsCache.update();
@@ -64,6 +65,206 @@ class GFS_LocalViewer
};
+//_________________________________________________________________________________________________________________________
+
+template<typename LFS>
+class BasisFN
+{
+public:
+ typedef typename LFS::Traits::FiniteElementType::Traits::LocalBasisType::Traits BTraits;
+ const std::vector<typename BTraits::RangeType>& phi;
+ const std::vector<typename BTraits::JacobianType>& gradphihat;
+ std::vector<typename BTraits::JacobianType> gradphi;
+
+ template<typename Pos, typename Cache, typename S>
+ BasisFN(const Cache& cache, const Pos& position, const LFS& femspace, S& s) :
+ phi(cache.evaluateFunction(position, femspace.finiteElement().localBasis())),
+ gradphihat(cache.evaluateJacobian(position, femspace.finiteElement().localBasis())),
+ gradphi(femspace.size())
+ {
+ for (size_t i=0; i<femspace.size(); i++)
+ {
+ s.mv(gradphihat[i][0],gradphi[i][0]);
+ }
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename CHLFS>
+class CHVariables
+{
+public:
+
+ const typename CHLFS::template Child<0>::Type& pfspace;
+ const typename CHLFS::template Child<1>::Type& xispace;
+ BasisFN<typename CHLFS::template Child<0>::Type> basis;
+
+ RF pf;
+ RF xi;
+ Dune::FieldVector<RF,dim> gradpf;
+ Dune::FieldVector<RF,dim> gradxi;
+
+ template<typename Pos, typename Cache, typename LV, typename S>
+ CHVariables(const Cache& cache, const Pos& position, const CHLFS& chLfs, const LV& lv, S& s) :
+ pfspace(child(chLfs,Dune::Indices::_0)),
+ xispace(child(chLfs,Dune::Indices::_1)),
+ basis(cache, position, pfspace, s),
+ pf(0.0),
+ xi(0.0),
+ gradpf(0.0),
+ gradxi(0.0)
+ {
+
+ for (size_t i=0; i<pfspace.size(); i++)
+ {
+ pf += lv[pfspace.localIndex(i)]*basis.phi[i];
+ xi += lv[xispace.localIndex(i)]*basis.phi[i];
+ gradpf.axpy(lv[pfspace.localIndex(i)],basis.gradphi[i][0]);
+ gradxi.axpy(lv[xispace.localIndex(i)],basis.gradphi[i][0]);
+ }
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename CHLFS>
+class CHVariables_U_PFOld_Extension : public CHVariables<RF,dim,CHLFS>
+{
+public:
+ const typename CHLFS::template Child<2>::Type& uspace;
+
+ RF pfOld;
+ Dune::FieldVector<RF,dim> gradpfOld;
+ RF u;
+ Dune::FieldVector<RF,dim> gradu;
+ RF uOld;
+
+
+ template<typename Pos, typename Cache, typename LV, typename S>
+ CHVariables_U_PFOld_Extension(const Cache& cache,
+ const Pos& position, const CHLFS& chLfs, const LV& lv, const LV& lvOld, S& s) :
+ CHVariables<RF,dim,CHLFS>(cache, position, chLfs, lv, s),
+ uspace(child(chLfs,Dune::Indices::_2)),
+ pfOld(0.0),
+ gradpfOld(0.0),
+ u(0.0),
+ gradu(0.0),
+ uOld(0.0)
+ {
+ for (size_t i=0; i<this->pfspace.size(); i++)
+ {
+ u += lv[uspace.localIndex(i)]*this->basis.phi[i];
+ gradu.axpy(lv[uspace.localIndex(i)],this->basis.gradphi[i][0]);
+ uOld += lvOld [uspace.localIndex(i)]*this->basis.phi[i];
+ pfOld += lvOld[this->pfspace.localIndex(i)]*this->basis.phi[i];
+ gradpfOld.axpy(lvOld[this->pfspace.localIndex(i)],this->basis.gradphi[i][0]);
+ }
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename NSLFS>
+class NSVariables
+{
+public:
+ const typename NSLFS::template Child<0>::Type& vspace;
+ const typename NSLFS::template Child<0>::Type::template Child<0>::Type& vfemspace;
+ BasisFN<typename NSLFS::template Child<0>::Type::template Child<0>::Type> vbasis;
+
+ Dune::FieldVector<RF,dim> v;
+ Dune::FieldMatrix<RF,dim,dim> jacv;
+ RF div_v;
+
+ template<typename Pos, typename vCache, typename LV, typename S>
+ NSVariables(const vCache& vcache, const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
+ vspace(child(nsLfs,Dune::Indices::_0)),
+ vfemspace(child(vspace,Dune::Indices::_0)),
+ vbasis(vcache, position, vfemspace, s),
+ v(0.0),
+ jacv(0.0),
+ div_v(0.0)
+ {
+ for(size_t k=0; k<dim; k++)
+ {
+ for (size_t i=0; i<vfemspace.size(); i++)
+ {
+ v[k] += lv[child(vspace,k).localIndex(i)]*vbasis.phi[i];
+ jacv[k].axpy(lv[child(vspace,k).localIndex(i)],vbasis.gradphi[i][0]);
+ }
+ div_v += jacv[k][k]; // div is trace of jacobian
+ }
+ }
+};
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename NSLFS>
+class NSVariables_P_Extension : public NSVariables<RF,dim,NSLFS>
+{
+public:
+ const typename NSLFS::template Child<1>::Type& pspace;
+ BasisFN<typename NSLFS::template Child<1>::Type> pbasis;
+
+ RF p;
+ Dune::FieldVector<RF,dim> gradp;
+
+ template<typename Pos, typename vCache, typename pCache, typename LV, typename S>
+ NSVariables_P_Extension(const vCache& vcache, const pCache& pcache,
+ const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
+ NSVariables<RF,dim,NSLFS>(vcache, position, nsLfs, lv, s),
+ pspace(child(nsLfs,Dune::Indices::_1)),
+ pbasis(pcache, position, pspace, s),
+ p(0.0),
+ gradp(0.0)
+ {
+ for (size_t i=0; i<pspace.size(); i++)
+ {
+ p += lv[pspace.localIndex(i)]*pbasis.phi[i];
+ gradp.axpy(lv[pspace.localIndex(i)],pbasis.gradphi[i][0]);
+ }
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename NSLFS>
+class NSVariables_P_VOld_Extension : public NSVariables_P_Extension<RF,dim,NSLFS>
+{
+public:
+ Dune::FieldVector<RF,dim> vOld;
+
+ template<typename Pos, typename vCache, typename pCache, typename LV, typename S>
+ NSVariables_P_VOld_Extension(const vCache& vcache, const pCache& pcache,
+ const Pos& position, const NSLFS& nsLfs, const LV& lv, const LV& lvOld, S& s) :
+ NSVariables_P_Extension<RF,dim,NSLFS>(vcache, pcache, position, nsLfs, lv, s),
+ vOld(0.0)
+ {
+
+ for (size_t i=0; i<this->vfemspace.size(); i++)
+ {
+ for(size_t k=0; k<dim; k++)
+ {
+ vOld[k] += lvOld[child(this->vspace,k).localIndex(i)]*this->vbasis.phi[i];
+ }
+ }
+ }
+};
+//_________________________________________________________________________________________________________________________
+
+
+template<typename CHGFS>
+using CHLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename CHGFS::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
+
+template<typename NSGFS>
+using VLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<0>::Type::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
+template<typename NSGFS>
+using PLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<1>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
// template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
// class FEM_CHNS_Base :
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
index 845a371ad4f66f9b2b7f7186d272ada7aed53f43..7ac9d0dcbc498e2a99c8fd6e4547de624edbc52a 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
@@ -8,59 +8,45 @@
#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_R_CahnHilliard :
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianBoundary<FEM_FS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyBoundary<FEM_FS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
{
private:
typedef typename CHContainer::field_type RF;
+ typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
- typedef typename CHGFS::template Child<0>::Type FEM;
- typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
- Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+ CHLocalBasisCache<CHGFS> chCache;
+ VLocalBasisCache<NSGFS> vCache;
- typedef typename NSGFS::template Child<0>::Type::template Child<0>::Type VFEM;
- typedef typename VFEM::Traits::FiniteElementType::Traits::LocalBasisType VLocalBasis;
- Dune::PDELab::LocalBasisCache<VLocalBasis> vcache;
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chOldLocal;
- Param& param; // parameter functions
+ typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+ mutable NSLView nsLocal;
- //FOR EVALUATING the old cahn-Hilliard values
- typedef Dune::PDELab::LocalFunctionSpace<CHGFS> CHLFS;
- mutable CHLFS chLfs;
- typedef Dune::PDELab::LFSIndexCache<CHLFS> CHLFSCache;
- mutable CHLFSCache chLfsCache;
- typedef typename CHContainer::template ConstLocalView<CHLFSCache> CHView;
- mutable CHView chView;
- mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chOldLocal;
-
- //FOR EVALUATING the old navier-Stokes values
- typedef Dune::PDELab::LocalFunctionSpace<NSGFS> NSLFS;
- mutable NSLFS nsLfs;
- typedef Dune::PDELab::LFSIndexCache<NSLFS> NSLFSCache;
- mutable NSLFSCache nsLfsCache;
- typedef typename NSContainer::template ConstLocalView<NSLFSCache> NSView;
- mutable NSView nsView;
- mutable std::vector<typename VLocalBasis::Traits::RangeFieldType> nsOldLocal;
-
- //FOR DETECTING ENTITY changes
- typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
- mutable ENT OldEnt;
+
+ Param& param; // parameter functions
public:
enum {doPatternVolume=true};
- //enum {doLambdaVolume=true};
enum {doAlphaVolume=true};
+ enum { doAlphaBoundary = true };
//constructor
- FEM_FS_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsVOld_) :
- param(param_),
- chLfs(chGfs_), chLfsCache(chLfs), chView(chVOld_), chOldLocal(chGfs_.maxLocalSize()),
- nsLfs(nsGfs_), nsLfsCache(nsLfs), nsView(nsVOld_), nsOldLocal(nsGfs_.maxLocalSize())
+ FEM_FS_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
+ chOldLocal(chGfs_,chVOld_),
+ nsLocal(nsGfs_,nsV),
+ param(param_)
{
}
@@ -68,150 +54,114 @@ public:
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
{
- //Left hand side, terms depending on the solution x
- using namespace Dune::Indices;
-
//Quadrature Rule
const int dim = EG::Entity::dimension;
const int order = 2; //TODO: Choose more educated
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
- //const auto& myordering = lfsv.gridFunctionSpace().ordering();
-
- const auto& pf1space = child(lfsv,_0);
- const auto& xi1space = child(lfsv,_1);
- const auto& uspace = child(lfsv,_2);
- const auto& femspace = child(lfsv,_0);
-
- const auto& vspace = child(nsLfs,_0);
- const auto& vfemspace = child(vspace,_0);
-
- if(eg.entity() != OldEnt)
- {
- OldEnt = eg.entity();
-
- chLfs.bind(eg.entity());
- chLfsCache.update();
- chView.bind(chLfsCache);
- chView.read(chOldLocal);
- chView.unbind();
-
- nsLfs.bind(eg.entity());
- nsLfsCache.update();
- nsView.bind(nsLfsCache);
-
- nsView.read_sub_container(vspace,nsOldLocal);
- nsView.unbind();
- }
+ chOldLocal.DoEvaluation(eg);
+ //chOldLocal.DoEvaluation(eg);
+ nsLocal.DoSubEvaluation(eg,child(nsLocal.lfs,Dune::Indices::_0));
for(const auto& ip : rule)
{
- // evaluate shape functions and gradient of shape functions
- auto& phi = cache.evaluateFunction(ip.position(), femspace.finiteElement().localBasis());
- auto& gradphihat = cache.evaluateJacobian(ip.position(), femspace.finiteElement().localBasis());
-
- auto& vphi = vcache.evaluateFunction(ip.position(), vfemspace.finiteElement().localBasis());
- auto& gradvphihat = vcache.evaluateJacobian(ip.position(), vfemspace.finiteElement().localBasis());
-
- // transform gradients of shape functions to real element
const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
- auto gradphi = makeJacobianContainer(femspace);
- auto gradvphi = makeJacobianContainer(vfemspace);
- for (size_t i=0; i<femspace.size(); i++)
- {
- S.mv(gradphihat[i][0],gradphi[i][0]);
- }
- for (size_t i=0; i<vfemspace.size(); i++)
- {
- S.mv(gradvphihat[i][0],gradvphi[i][0]);
- }
- //compute current values of the unknowns
- Dune::FieldVector<RF,dim> v(0.0);
- Dune::FieldMatrix<RF,dim,dim> jacv(0.0);
-
- RF pf1=0.0;
- RF pf1Old=0.0;
- RF xi1=0.0;
-
- RF u=0.0;
- RF uOld=0.0;
-
- Dune::FieldVector<RF,dim> gradpf1(0.0);
- Dune::FieldVector<RF,dim> gradpf1Old(0.0);
- Dune::FieldVector<RF,dim> gradxi1(0.0);
-
- Dune::FieldVector<RF,dim> gradu(0.0);
-
- for (size_t i=0; i<femspace.size(); i++)
- {
-
- gradpf1.axpy(x(pf1space,i),gradphi[i][0]);
- pf1 += x(pf1space,i)*phi[i];
- gradxi1.axpy(x(xi1space,i),gradphi[i][0]);
- xi1 += x(xi1space,i)*phi[i];
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
- gradu.axpy(x(uspace,i),gradphi[i][0]);
- u += x(uspace,i)*phi[i];
+ 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);
- pf1Old += chOldLocal[pf1space.localIndex(i)]*phi[i];
- gradpf1Old.axpy(chOldLocal[pf1space.localIndex(i)],gradphi[i][0]);
- uOld += chOldLocal[uspace.localIndex(i)]*phi[i];
- }
- for (size_t i=0; i<vfemspace.size(); i++)
- {
- for(size_t k=0; k<dim; k++)
- {
- v[k] += nsOldLocal[child(vspace,k).localIndex(i)]*vphi[i];
- jacv[k].axpy(nsOldLocal[child(vspace,k).localIndex(i)],gradvphi[i][0]);
- }
- }
-
- RF div_v = 0;
- for(int k=0; k<dim; ++k)
- {
- div_v += jacv[k][k]; // divergence of v is the trace of the velocity jacobian
- }
-
- RF pf_c = pf1;
- RF pf_c_Old = pf1Old;
+ RF pf_c = ch.pf;
+ RF pf_c_Old = ch.pfOld;
RF pf_c_reg = pf_c+param.delta;
RF pf_c_reg_Old = pf_c_Old + param.delta;
- auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
- for (size_t i=0; i<femspace.size(); i++)
+ for (size_t i=0; i<ch.pfspace.size(); i++)
{
- RF R1 = -param.dw.eval_q(pf1Old)/param.eps*
- (param.reaction.eval(u)+param.phys.curvatureAlpha*xi1)*phi[i]; //TODO Scale with sigma13
- RF J1 = param.mobility.eval(pf1)*param.eps * (gradxi1*gradphi[i][0]);
+ RF R1 = -std::max(param.dw.eval_q(ch.pf) - 0.05, 0.0)/param.eps*
+ (param.reaction.eval(ch.u)+param.phys.curvatureAlpha*ch.xi)*ch.basis.phi[i]; //TODO Scale with sigma13
+ RF J1 = param.mobility.eval(ch.pf)*param.eps * (ch.gradxi*ch.basis.gradphi[i][0]);
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
- r.accumulate(pf1space,i,factor*(
- (pf1-pf1Old)*phi[i]/param.time.dt
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+J1
- //+(pf1*div_v + (gradpf1*v))*phi[i]
+ //+(pf1*div_v + (gradpf1*v))*phi[i] TODO ADD THIS AGAIN WHEN BC ARE FIXED
-R1
));
- r.accumulate(xi1space,i,factor*(
- xi1*phi[i]
- - 1/param.eps* param.dw.eval_d(pf1) *phi[i]
- - param.eps * (gradpf1*gradphi[i][0])
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.dw.eval_d(ch.pf) *ch.basis.phi[i]
+ - param.eps * (ch.gradpf*ch.basis.gradphi[i][0])
));
// 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(uspace,i,factor*(
+ r.accumulate(ch.uspace,i,factor*(
//(u-param.uAst)*phi[i]/param.dt //For Testing
- (pf_c_reg-pf_c_reg_Old)*(u-param.phys.uAst)*phi[i]/param.time.dt
- + pf_c_reg_Old*(u-uOld)*phi[i]/param.time.dt
- + param.phys.D*pf_c_reg_Old*(gradu*gradphi[i][0])
- + J1*(u-param.phys.uAst)
- + (pf_c*(u-param.phys.uAst)*div_v + (u-param.phys.uAst)*(gradpf1*v) + pf_c*(gradu*v) )*phi[i]
+ (pf_c_reg-pf_c_reg_Old)*(ch.u-param.phys.uAst)*ch.basis.phi[i]/param.time.dt
+ + 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])
+ + J1*(ch.u-param.phys.uAst)
+ + (pf_c*(ch.u-param.phys.uAst)*ns.div_v + (ch.u-param.phys.uAst)*(ch.gradpf*ns.v) + pf_c*(ch.gradu*ns.v) )*ch.basis.phi[i]
));
}
}
}
+ // boundary integral
+ template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_boundary (const IG& ig,
+ const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
+ R& r_s) const
+ {
+ const int dim=2;
+ // Define types
+ using size_type = typename LFSV::Traits::SizeType;
+
+ const auto& pfspace = child(lfsv_s,Dune::Indices::_0);
+
+ // Reference to the inside cell
+ const auto& cell_inside = ig.inside();
+
+ // Get geometry
+ auto geo = ig.geometry();
+ auto ref_el_ig = referenceElement(geo);
+ auto face_center_ig = ref_el_ig.position(0,0);
+
+ // Get geometry of intersection in local coordinates of cell_inside
+ auto geo_in_inside = ig.geometryInInside();
+
+ auto geo_inside = cell_inside.geometry();
+
+ // evaluate boundary condition type
+ auto ref_el = referenceElement(geo_in_inside);
+ auto local_face_center = ref_el.position(0,0);
+ auto intersection = ig.intersection();
+
+ Dune::FieldVector<RF, dim> xg = ig.geometry().global( face_center_ig);
+
+ if(xg[0]>1e-7) return; //Only non-dirichlet at left brdy
+
+ // loop over quadrature points and integrate normal flux
+ auto intorder = 4;
+ for (const auto& ip : quadratureRule(geo,intorder))
+ {
+ auto local = geo_in_inside.global(ip.position());
+ // position of quadrature point in local coordinates of element
+ const auto& S = geo_inside.jacobianInverseTransposed(local);
+
+ auto n = ig.unitOuterNormal(ip.position());
+
+ CHVariables<RF,dim,LFSV> ch(chCache, local, lfsv_s, x_s.base(), S);
+
+ auto factor = ip.weight()*geo.integrationElement(ip.position());
+ for (size_type i=0; i<ch.pfspace.size(); i++)
+ r_s.accumulate(ch.xispace,i,0);//-(ch.gradpf*n)*ch.basis.phi[i]*factor);
+
+ }
+ }
};
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
index fe559dabd1f2102c43340ced2ed8c7bcb823f5e3..d4551beb76f90c383eec4aeeb8b84207b09fd0ea 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
@@ -22,44 +22,21 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
{
private:
typedef typename CHContainer::field_type RF;
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
- typedef typename CHGFS::template Child<0>::Type FEM;
- typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
- Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+ CHLocalBasisCache<CHGFS> chCache;
+ VLocalBasisCache<NSGFS> vCache;
+ PLocalBasisCache<NSGFS> pCache;
- typedef typename NSGFS::template Child<0>::Type::template Child<0>::Type VFEM;
- typedef typename VFEM::Traits::FiniteElementType::Traits::LocalBasisType VLocalBasis;
- Dune::PDELab::LocalBasisCache<VLocalBasis> vcache;
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chLocal;
+ //mutable CHLView chOldLocal;
- typedef typename NSGFS::template Child<1>::Type PFEM;
- typedef typename PFEM::Traits::FiniteElementType::Traits::LocalBasisType PLocalBasis;
- Dune::PDELab::LocalBasisCache<PLocalBasis> pcache;
+ typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+ mutable NSLView nsOldLocal;
Param& param; // parameter functions
- //FOR EVALUATING the new and old cahn-Hilliard values
- typedef Dune::PDELab::LocalFunctionSpace<CHGFS> CHLFS;
- mutable CHLFS chLfs;
- typedef Dune::PDELab::LFSIndexCache<CHLFS> CHLFSCache;
- mutable CHLFSCache chLfsCache;
- typedef typename CHContainer::template ConstLocalView<CHLFSCache> CHView;
- mutable CHView chView;
- mutable CHView chOldView;
- mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chLocal;
- mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chOldLocal;
-
- //FOR EVALUATING the old navier-Stokes values
- typedef Dune::PDELab::LocalFunctionSpace<NSGFS> NSLFS;
- mutable NSLFS nsLfs;
- typedef Dune::PDELab::LFSIndexCache<NSLFS> NSLFSCache;
- mutable NSLFSCache nsLfsCache;
- typedef typename NSContainer::template ConstLocalView<NSLFSCache> NSView;
- mutable NSView nsView;
- mutable std::vector<typename VLocalBasis::Traits::RangeFieldType> nsOldLocal;
-
- //FOR DETECTING ENTITY changes
- typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
- mutable ENT OldEnt;
public:
//! Boundary condition indicator type
@@ -79,12 +56,12 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
//TODO enum { doLambdaBoundary = true };
- FEM_FS_CHNS_R_InstatStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
+ FEM_FS_CHNS_R_InstatStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
NSGFS& nsGfs_, NSContainer& nsVOld_) :
- param(param_),
- chLfs(chGfs_), chLfsCache(chLfs), chView(chV_), chOldView(chVOld_),
- chLocal(chGfs_.maxLocalSize()), chOldLocal(chGfs_.maxLocalSize()),
- nsLfs(nsGfs_), nsLfsCache(nsLfs), nsView(nsVOld_), nsOldLocal(nsGfs_.maxLocalSize())
+ chLocal(chGfs_,chV_),
+ //chOldLocal(chGfs_,chVOld_),
+ nsOldLocal(nsGfs_,nsVOld_),
+ param(param_)
{
}
@@ -92,152 +69,55 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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
{
- //Left hand side, terms depending on the solution x
- using namespace Dune::Indices;
-
//Quadrature Rule
const int dim = EG::Entity::dimension;
const int order = 4; //TODO: Choose more educated
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
- //const auto& myordering = lfsv.gridFunctionSpace().ordering();
-
- const auto& pfspace = child(chLfs,_0);
- const auto& xispace = child(chLfs,_1);
- const auto& femspace = child(chLfs,_0);
-
- const auto& vspace = child(lfsv,_0);
- const auto& pspace = child(lfsv,_1);
- const auto& vfemspace = child(vspace,_0);
-
- if(eg.entity() != OldEnt)
- {
- OldEnt = eg.entity();
- chLfs.bind(eg.entity());
- chLfsCache.update();
-
- chView.bind(chLfsCache);
- chView.read(chLocal);
- chView.unbind();
-
- nsLfs.bind(eg.entity());
- nsLfsCache.update();
-
- nsView.bind(nsLfsCache);
- nsView.read_sub_container(vspace,nsOldLocal);
- nsView.unbind();
- }
+ chLocal.DoEvaluation(eg);
+ //chOldLocal.DoEvaluation(eg);
+ nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
for(const auto& ip : rule)
{
- // evaluate shape functions and gradient of shape functions
- auto& phi = cache.evaluateFunction(ip.position(), femspace.finiteElement().localBasis());
- auto& gradphihat = cache.evaluateJacobian(ip.position(), femspace.finiteElement().localBasis());
-
- auto& vphi = vcache.evaluateFunction(ip.position(), vfemspace.finiteElement().localBasis());
- auto& gradvphihat = vcache.evaluateJacobian(ip.position(), vfemspace.finiteElement().localBasis());
-
- auto& pphi = pcache.evaluateFunction(ip.position(), pspace.finiteElement().localBasis());
- auto& gradpphihat = pcache.evaluateJacobian(ip.position(), pspace.finiteElement().localBasis());
-
- // transform gradients of shape functions to real element
const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
- auto gradphi = makeJacobianContainer(femspace);
- auto gradvphi = makeJacobianContainer(vfemspace);
- auto gradpphi = makeJacobianContainer(pspace);
- for (size_t i=0; i<femspace.size(); i++)
- {
- S.mv(gradphihat[i][0],gradphi[i][0]);
- }
- for (size_t i=0; i<vfemspace.size(); i++)
- {
- S.mv(gradvphihat[i][0],gradvphi[i][0]);
- }
- for (size_t i=0; i<pspace.size(); i++)
- {
- S.mv(gradpphihat[i][0],gradpphi[i][0]);
- }
-
- //compute current values of the unknowns
- Dune::FieldVector<RF,dim> v(0.0);
- Dune::FieldVector<RF,dim> vOld(0.0);
- Dune::FieldMatrix<RF,dim,dim> jacv(0.0);
- RF p=0.0;
- Dune::FieldVector<RF,dim> gradp(0.0);
- RF pf=0.0;
- Dune::FieldVector<RF,dim> gradpf(0.0);
- Dune::FieldVector<RF,dim> gradxi(0.0);
-
- for (size_t i=0; i<pspace.size(); i++)
- {
- p += x(pspace,i)*pphi[i];
- gradp.axpy(chLocal[pspace.localIndex(i)],gradpphi[i][0]);
- }
-
- for (size_t i=0; i<femspace.size(); i++)
- {
- pf += chLocal[pfspace.localIndex(i)]*phi[i];
- gradpf.axpy(chLocal[pfspace.localIndex(i)],gradphi[i][0]);
- gradxi.axpy(chLocal[xispace.localIndex(i)],gradphi[i][0]);
- }
- for (size_t i=0; i<vfemspace.size(); i++)
- {
- for(size_t k=0; k<dim; k++)
- {
- v[k] += x(child(vspace,k),i)*vphi[i];
- vOld[k] += nsOldLocal[child(vspace,k).localIndex(i)]*vphi[i];
- jacv[k].axpy(x(child(vspace,k),i),gradvphi[i][0]);
- }
- }
-
- Dune::FieldVector<RF,dim> g(0.0);
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
- //Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
- //g[1]= 20*xg[0]*(1-xg[0]);
- g[1] = -5;
+ CHVariables<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+ NSVariables_P_VOld_Extension<RF,dim,LFSV> ns(vCache, pCache, ip.position(), lfsv, x.base(), nsOldLocal.lv, S);
- RF rho_f=param.phys.rho_1*pf;
+ RF rho_f=param.phys.rho_1*ch.pf;
RF rho_f_reg= rho_f + param.delta;
- RF pf_f = pf + 2*param.delta*(1-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),gradpf);
+ gradpf_f.axpy((1-2*param.delta),ch.gradpf);
RF mu_reg = pf_f*param.phys.mu;
- RF div_v = 0;
- for(int k=0; k<dim; ++k)
- {
- div_v += jacv[k][k]; // divergence of v is the trace of the velocity jacobian
- }
-
- auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
- for (size_t i=0; i<pspace.size(); i++)
+ for (size_t i=0; i<ns.pspace.size(); i++)
{
//r.accumulate(pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);
- r.accumulate(pspace,i, -1.0 * (pf_f*div_v + (gradpf_f*v) )* pphi[i] * factor);
+ r.accumulate(ns.pspace,i, -1.0 * (pf_f*ns.div_v + (gradpf_f*ns.v) )* ns.pbasis.phi[i] * factor);
}
- for (size_t i=0; i<vfemspace.size(); i++)
+
+ for (size_t i=0; i<ns.vfemspace.size(); i++)
{
for(int k=0; k<dim; k++)
{
//SIMPLIFY S=0, TODO: Implement this correctly
-
- const auto v_nabla_v = v * jacv[k];
-
-
- r.accumulate(child(vspace,k),i, ( rho_f_reg*(v[k]-vOld[k])*vphi[i]/param.time.dt //TimeEvolution
- - p * (pf_f*gradvphi[i][0][k] + vphi[i]*gradpf_f[k]) //Pressure in weak form
+ r.accumulate(child(ns.vspace,k),i, ( rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt //TimeEvolution
+ - ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k]) //Pressure in weak form
//+ pf_f * (gradp[k] * vphi[i]) //Pressure in strong form
- + mu_reg * (jacv[k]*gradvphi[i][0]) //Viscosity
- + param.vDis.eval(pf_f)*v[k]*vphi[i] //Velocity Dissipation in Solid
+ + mu_reg * (ns.jacv[k]*ns.vbasis.gradphi[i][0]) //Viscosity
+ + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i] //Velocity Dissipation in Solid
//- rho*g[k]*vphi[i] //Gravitation
//+ S //Surface Tension
//+ Reaction TODO
)*factor);
// integrate (grad u)^T * grad phi_i //Transpose for Viscosity
for(int dd=0; dd<dim; ++dd)
- r.accumulate(child(vspace,k),i, param.phys.mu * (jacv[dd][k] * (gradvphi[i][0][dd])) * factor);
+ r.accumulate(child(ns.vspace,k),i, param.phys.mu * (ns.jacv[dd][k] * (ns.vbasis.gradphi[i][0][dd])) * factor);
}
}
}
diff --git a/src/test_2pchns_r/2p_driver.hh b/src/test_2pchns_r/2p_driver.hh
index f949feb8c7297f5421c930b251c4a26ba3777350..c636cd8d830765e4599fa7e0e196ed7af4a374c9 100644
--- a/src/test_2pchns_r/2p_driver.hh
+++ b/src/test_2pchns_r/2p_driver.hh
@@ -37,7 +37,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
typedef Params_fs_r< RF,
- DoubleWell_poly<RF>,
+ DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,//DoubleWell_poly<RF>,
Mobility_Const<RF>,
Reaction_Quadratic_Limited<RF>,
VelDissipation_Quadratic<RF> > Parameters;
@@ -58,7 +58,8 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
//_________________________________________________________________________________________________________________________
//Constraints
- typedef Bdry_fs_chns_r<dim, InflowOutflow_v,
+ typedef Bdry_fs_chns_r_compositeV_2d<InflowOutflow_v,
+ Dune::PDELab::DirichletConstraintsParameters,
InflowOutflow_p,
InflowOutflow_phi,
Dune::PDELab::NoDirichletConstraintsParameters,
@@ -79,7 +80,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
typedef Ini_fs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
ZeroInitial<GV,RF,Parameters>,
- PhiInitial_Sphere<GV,RF,Parameters>,
+ PhiInitial_Rectangle<GV,RF,Parameters>,
ZeroInitial<GV,RF,Parameters>,
UInitial_Inflow<GV,RF,Parameters> > Initial;
Initial initial(gv,param);
@@ -89,7 +90,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
//____________________________________________________________________________________________________________________
// Refine Initial Data
- for(int i=0; i<4; i++)
+ for(int i=0; i<ptree.get("Refinement.maxRefineLevel",2); i++)
{
// mark elements for refinement
typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
@@ -117,7 +118,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
// LocalOperators
typedef FEM_FS_CHNS_R_InstatStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
- NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+ NS_LOP nsLop(param, gfs.ch, chV, gfs.ns, nsVOld);
typedef FEM_FS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
diff --git a/src/test_2pchns_r/2pchns_r.ini b/src/test_2pchns_r/2pchns_r.ini
index fcdb12312ce985341d944de91a1d63ba7744131d..c9bb158fb539ad75c78c2e0895ea22239fa10a76 100644
--- a/src/test_2pchns_r/2pchns_r.ini
+++ b/src/test_2pchns_r/2pchns_r.ini
@@ -1,27 +1,29 @@
-[physics]
+[Physics]
uAst = 2
-D = 0.005
+D = 0.01
rho1 = 1.0
rho2 = 1.0
mu=0.01
-ReactionRate=0 #=3
+ReactionRate=0.1
#ReactionLimiter = 1000
-#VelDissipationRate = 10
-CurvatureAlpha=0.03
+VelDissipationRate = 100
+CurvatureAlpha=0.001
[domain]
level = 0
#filename = grids/unitsquare.msh
#filename = grids/square_27k_del.msh
#filename = grids/square_7k_del.msh
-filename = grids/rectangle_3k_per.msh
+#filename = grids/rectangle_3k_per.msh
#filename = grids/rectangle_53k_per.msh
-#filename = grids/square_2k_del_per.msh
+filename = grids/square_2k_del_per.msh
[Phasefield]
eps=0.005
-delta=0.02
+delta=0.01
Mobility=0.005
+DoubleWellDelta=0.005
+DoubleWellGamma=0.0045
[Time]
tMax = 30