diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 86f1a1d6cc0dbe3b3c9a677e8ae05062c3d8b0cf..6d8b0a0bc1708be66da1051234eadf2361b6e2cf 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -232,8 +232,7 @@ public:
Number eval_mu (const Number& pf_f) const
{
- //return pf_f*mu+(1-pf_f)*muS;
- return 1/(pf_f*1/mu + (1-pf_f)*1/muS); //TODO Evaluate differences
+ return 1/(pf_f*1/mu + (1-pf_f)*1/muS);
}
};
@@ -396,87 +395,3 @@ public:
{
}
};
-//
-// template<typename Number>
-// class Problem
-// {
-// Number t;
-//
-//
-// public:
-// typedef Number value_type;
-// Number eps, delta;
-//
-// Number dt;
-//
-// Number u_regularization;
-//
-//
-// Number M;
-// Number M_reduced;
-// Number fscaling;
-// Number SurfaceTensionScaling;
-//
-// //Number rho (const Number& x) {return rho_0;}
-//
-//
-//
-//
-//
-// // Constructor takes parameter
-// Problem (const Number& dt_, const Dune::ParameterTree& ptree)
-// : t(0.0)
-// {
-// eps = ptree.get("problem.eps",(Number)0.1);
-// delta = ptree.get("problem.delta",(Number)0.1);
-// u_regularization = ptree.get("problem.URegularization",(Number)0.01);
-//
-//
-// fscaling = ptree.get("phsyics.fscaling",(Number)1);
-//
-//
-// dt=dt_;
-// }
-//
-
-//
-//
-// Number f (const Number& u) const
-// {
-// if(u<0)
-// return fscaling*(-1.0+1000*u);
-// else
-// return fscaling*(u*u-1);
-// }
-//
-// Number q (const Number& pf) const
-// {
-// if (pf<0 || pf>1)
-// return 0;
-// return pf*(1-pf);
-// }
-//
-// Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
-// {
-// if (pf1<0 || pf3<0)
-// return 0;
-// return pf1*pf3;
-// }
-//
-// Number Diffusion(const Number& pf_c)
-// {
-// return D*fmax(1e-3,(pf_c-0.05)/0.95);
-// }
-//
-// Number VelDissipation(const Number& pf_f)
-// {
-// return 10*(1-pf_f)*(1-pf_f);
-// }
-//
-// //! Set time in instationary case
-// void setTime (Number t_)
-// {
-// t = t_;
-// }
-//
-// };
diff --git a/dune/phasefield/base_potentials.hh b/dune/phasefield/base_potentials.hh
index dd199383184392109f798e20d70a8f0a4647ed20..b5f70dc41a95521a45b7549e2568cc6092a4c953 100644
--- a/dune/phasefield/base_potentials.hh
+++ b/dune/phasefield/base_potentials.hh
@@ -209,7 +209,9 @@ public:
//______________________________________________________________________________________________________
template<typename Number>
-class Limiter_FakeDiverging
+class Limiter_FakeDiverging //A potential that uses the limiter delta x^2/(delta+x) up to the point gamma
+ //and extends the limiter by a linear function from there on
+ //use with 0 < delta - gamma << 1 for better convergence with newton than Limiter_Diverging
{
public:
const Number delta; //Where the potential would be diverging
diff --git a/dune/phasefield/ff_chns.hh b/dune/phasefield/ff_chns.hh
index a8414da7378bc4e3d41f713a94b0454176b77e9e..3a846df22a4f0e1ef2a4e372d6cc3e914222c82f 100644
--- a/dune/phasefield/ff_chns.hh
+++ b/dune/phasefield/ff_chns.hh
@@ -179,3 +179,34 @@ public:
{
}
};
+
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu>
+class Bdry_ff_chns_r_compositeV_2d_Explicit_Domain
+{
+public:
+ BdryV1 bdryV1;
+ BdryV2 bdryV2;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ template<typename DomainSize>
+ Bdry_ff_chns_r_compositeV_2d_Explicit_Domain(DomainSize ds) :
+ bdryV1(ds),
+ bdryV2(ds),
+ bdryP(ds),
+ bdryPhi(ds),
+ bdryMu(ds),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu)
+ {
+ }
+};
diff --git a/dune/phasefield/fs_chns_r.hh b/dune/phasefield/fs_chns_r.hh
index e2e04a39f17544c9bfed72addc0e9dcee5668d66..a53113841d5d6cef6fd81a8b1562467346a5f014 100644
--- a/dune/phasefield/fs_chns_r.hh
+++ b/dune/phasefield/fs_chns_r.hh
@@ -229,3 +229,36 @@ public:
{
}
};
+
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryU>
+class Bdry_fs_chns_r_compositeV_2d_Explicit_Domain
+{
+public:
+ 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;
+
+ template<typename DomainSize>
+ Bdry_fs_chns_r_compositeV_2d_Explicit_Domain(DomainSize ds) :
+ bdryV1(ds),
+ bdryV2(ds),
+ bdryP(ds),
+ bdryPhi(ds),
+ bdryMu(ds),
+ bdryU(ds),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryU)
+ {
+ }
+};
diff --git a/dune/phasefield/localoperator/fem_2p_cahnhilliard.hh b/dune/phasefield/localoperator/fem_2p_cahnhilliard.hh
index e7745a5c04f5094d604237fe49e9dc4a5ebcba55..de57851cfd761a8e92e4bd5effd7363ce0970ba5 100644
--- a/dune/phasefield/localoperator/fem_2p_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_2p_cahnhilliard.hh
@@ -47,7 +47,7 @@ public:
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
+ const int order = 5;
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
diff --git a/dune/phasefield/localoperator/fem_3p_cahnhilliard.hh b/dune/phasefield/localoperator/fem_3p_cahnhilliard.hh
index 98b603b26d9d1ced0230e6fc4eba4e6077a5694d..0225ecc5e74f950b540af2428e3d2f0941166587 100644
--- a/dune/phasefield/localoperator/fem_3p_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_3p_cahnhilliard.hh
@@ -11,7 +11,6 @@
#include<dune/phasefield/localoperator/fem_base_chns.hh>
-//Parameter class (for functions), FEM space and RangeField
template<typename Param, typename CHGFS, typename CHContainer>
class FEM_3P_CahnHilliard :
public Dune::PDELab::FullVolumePattern,
@@ -45,9 +44,8 @@ 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
{
- //Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
+ const int order = 5; //Quadrature Rule
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
diff --git a/dune/phasefield/localoperator/fem_base_chns.hh b/dune/phasefield/localoperator/fem_base_chns.hh
index 98de02b985b31252795c331370643fac05ed6bc2..41eaf878b3b4cca7639cfdc5874e36bf849d940f 100644
--- a/dune/phasefield/localoperator/fem_base_chns.hh
+++ b/dune/phasefield/localoperator/fem_base_chns.hh
@@ -493,6 +493,8 @@ public:
//_________________________________________________________________________________________________________________________
+// Local Basis Caches used for Operator Splitting
+
template<typename CHGFS>
using CHLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename CHGFS::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
@@ -506,11 +508,13 @@ template<typename NSGFS>
using PLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<1>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
+// Local Basis Caches used for Monolithic Structure
+
template<typename CHNSGFS>
using CH_CHNSLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<0>::Type::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
-
template<typename CHNSGFS>
using V_CHNSLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<1>::Type::template Child<0>::Type::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
@@ -518,235 +522,3 @@ using V_CHNSLocalBasisCache =
template<typename CHNSGFS>
using P_CHNSLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<1>::Type::template Child<1>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
-
-// template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
-// class FEM_CHNS_Base :
-// public Dune::PDELab::FullVolumePattern,
-// public Dune::PDELab::LocalOperatorDefaultFlags,
-// public Dune::PDELab::NumericalJacobianVolume<FEM_CHNS_Base<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
-// public Dune::PDELab::NumericalJacobianApplyVolume<FEM_CHNS_Base<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
-// {
-// private:
-// typedef typename CHContainer::field_type RF;
-//
-// typedef typename CHGFS::template Child<0>::Type FEM;
-// typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
-// Dune::PDELab::LocalBasisCache<LocalBasis> cache;
-//
-// 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 typename NSGFS::template Child<1>::Type PFEM;
-// typedef typename PFEM::Traits::FiniteElementType::Traits::LocalBasisType PLocalBasis;
-// Dune::PDELab::LocalBasisCache<PLocalBasis> pcache;
-//
-// 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
-// // using BC = StokesBoundaryCondition; // ??
-//
-//
-// // static const bool navier = P::assemble_navier;
-// // static const bool full_tensor = P::assemble_full_tensor;
-//
-// // pattern assembly flags
-// enum { doPatternVolume = true };
-//
-// // residual assembly flags
-// enum { doAlphaVolume = true };
-// //enum { doLambdaVolume = true };
-//
-// //TODO enum { doLambdaBoundary = true };
-//
-//
-// FEM_CHNS_Base (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
-// 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())
-// {
-// }
-//
-//
-// 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();
-// }
-//
-// 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);
-//
-// //Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
-// //g[1]= 20*xg[0]*(1-xg[0]);
-// g[1] = -5;
-//
-// RF rho_f=param.phys.rho_1*pf;
-// RF rho_f_reg= rho_f + param.delta;
-//
-// RF pf_f = pf + 2*param.delta*(1-pf);
-// Dune::FieldVector<RF,dim> gradpf_f(0.0);
-// gradpf_f.axpy((1-2*param.delta),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++)
-// {
-// //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);
-// }
-// for (size_t i=0; i<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
-// //+ 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
-// //- 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);
-// }
-// }
-// }
-// }
-//
-// };
diff --git a/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
index ab36cc779f775be3a184a63b0eccb635f3169a07..b1023899f2a20626fa6c550f6d914efe422adea2 100644
--- a/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
@@ -8,7 +8,6 @@
#include<dune/pdelab/localoperator/variablefactories.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
-//Parameter class (for functions), FEM space and RangeField
template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
class FEM_FF_CHNS_CahnHilliard :
public Dune::PDELab::FullVolumePattern,
@@ -73,9 +72,8 @@ public:
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 2; //TODO: Choose more educated
+ const int order = 2; //Integration order impacts performance & error
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
- //const auto& myordering = lfsv.gridFunctionSpace().ordering();
const auto& pfspace = child(lfsv,_0);
const auto& xispace = child(lfsv,_1);
diff --git a/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh b/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
index 9461fda8c67e496e21633383b43ed40f58ffc659..bc69754a0475154ad1e0d23c2231ac65a9f5d19d 100644
--- a/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
+++ b/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
@@ -2,7 +2,6 @@
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/finiteelement/localbasiscache.hh>
-//#include<dune/pdelab/localoperator/idefault.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
@@ -56,21 +55,11 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
mutable ENT OldEnt;
public:
- //! Boundary condition indicator type
- // using BC = StokesBoundaryCondition; // ??
-
-
-// static const bool navier = P::assemble_navier;
-// static const bool full_tensor = P::assemble_full_tensor;
-
- // pattern assembly flags
+ // pattern assembly flags
enum { doPatternVolume = true };
// residual assembly flags
enum { doAlphaVolume = true };
- //enum { doLambdaVolume = true };
-
- //TODO enum { doLambdaBoundary = true };
FEM_FF_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
index 834972d7f0fc3e23da52593d11bf9260236dc09d..e44c367e42884d11c30e2bf4a4dbd79d7813d116 100644
--- a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
@@ -82,12 +82,12 @@ public:
RF R1 = -Rspeed/param.eps*
(param.reaction.eval(ch.u)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i];
RF J1 = param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
-
+
// 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
+J1
- //+(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
+ //+(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i] //strong transport
-ch.pf*(ns.v*ch.basis.gradphi[i][0]) //weak transport
-R1
));
@@ -95,7 +95,7 @@ public:
r.accumulate(ch.pf2space,i,factor*(
(ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0])
- //+(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
+ //+(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i] //strong transport
-ch.pf2*(ns.v*ch.basis.gradphi[i][0])
));
@@ -116,20 +116,20 @@ public:
Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
r.accumulate(ch.uspace,i,factor*(
- //(u-param.uAst)*phi[i]/param.dt //For Testing
-
(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*(u-param.phys.uAst)*(v*gradphi[i][0])
- + (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]
+ //- pf_c*(u-param.phys.uAst)*(v*gradphi[i][0]) //Weak form of transport
+ + (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] //Strong form of transport
));
- if(xg[0]<=0.1)
+
+ //Source term in the domain
+ /*if(xg[0]<=0.1)
{
r.accumulate(ch.uspace,i,factor*(
10*(ch.u-1.5)*ch.basis.phi[i] ));
- }
+ }*/
}
}
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh
index a8c06c5fa0d62d03b8778fca76c9a23cf492884b..ed13e35fb3698156095f31276882046742bf596b 100644
--- a/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh
@@ -69,11 +69,10 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
+ const int order = 5;
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chLocal.DoEvaluation(eg);
- //chOldLocal.DoEvaluation(eg);
nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
for(const auto& ip : rule)
@@ -108,14 +107,9 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
{
// integrate: ( (rho+rhoOld)/2 * v_n+1 - rhoOld * v_n)/dt * phi + nu * Dv_n+1 : D phi - p_n+1 div phi
- //std::cout << (param.nu * (jacv[k]*gradvphi[i][0]) - p * gradvphi[i][0][k]
- //- param.rho*g[k]*vphi[i] )*factor << std::endl;
-
- //RF lambda = 3;
- //const auto v_nabla_v = v * jacv[k];
RF gradxi3k = -param.tw.Sigma3*(ch.gradxi[k]/param.tw.Sigma1 + ch.gradxi2[k]/param.tw.Sigma2);
RF S = 0;
- S += - ch.xi2*(ch.gradpf[k]*pf_f - ch.pf * gradpf_f[k])/(pf_f*pf_f); //TODO!!!!!!! Other pf_f scaling!!!!!
+ S += - ch.xi2*(ch.gradpf[k]*pf_f - ch.pf * gradpf_f[k])/(pf_f*pf_f);
S += - ch.xi *(ch.gradpf2[k]*pf_f - ch.pf2 * gradpf_f[k])/(pf_f*pf_f);
S += - 2*param.delta*pf3*(gradxi3k-ch.gradxi2[k]-ch.gradxi[k]);
diff --git a/dune/phasefield/localoperator/fem_fs_chns_driven_thinplateflow.hh b/dune/phasefield/localoperator/fem_fs_chns_driven_thinplateflow.hh
deleted file mode 100644
index ef04bc23e3351dda23097398a1d1996964daeba9..0000000000000000000000000000000000000000
--- a/dune/phasefield/localoperator/fem_fs_chns_driven_thinplateflow.hh
+++ /dev/null
@@ -1,130 +0,0 @@
-#pragma once
-
-#include<dune/pdelab/common/quadraturerules.hh>
-#include<dune/pdelab/finiteelement/localbasiscache.hh>
-//#include<dune/pdelab/localoperator/idefault.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>
-
-
-
-template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
- class FEM_FS_CHNS_Driven_ThinPlateFlow :
- public Dune::PDELab::FullVolumePattern,
- public Dune::PDELab::LocalOperatorDefaultFlags,
- public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_Driven_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_Driven_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
- {
- private:
- typedef typename CHContainer::field_type RF;
- typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
-
- CHLocalBasisCache<CHGFS> chCache;
- VLocalBasisCache<NSGFS> vCache;
- PLocalBasisCache<NSGFS> pCache;
-
- typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
- mutable CHLView chLocal;
- //mutable CHLView chOldLocal;
-
- typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
- mutable NSLView nsOldLocal;
-
- Param& param; // parameter functions
-
-
- public:
- //! Boundary condition indicator type
- // using BC = StokesBoundaryCondition; // ??
-
-
-// static const bool navier = P::assemble_navier;
-// static const bool full_tensor = P::assemble_full_tensor;
-
- // pattern assembly flags
- enum { doPatternVolume = true };
-
- // residual assembly flags
- enum { doAlphaVolume = true };
- //enum { doLambdaVolume = true };
-
- //TODO enum { doLambdaBoundary = true };
-
-
- FEM_FS_CHNS_Driven_ThinPlateFlow (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
- NSGFS& nsGfs_, NSContainer& nsVOld_) :
- chLocal(chGfs_,chV_),
- //chOldLocal(chGfs_,chVOld_),
- nsOldLocal(nsGfs_,nsVOld_),
- 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);
-
- chLocal.DoEvaluation(eg);
- //chOldLocal.DoEvaluation(eg);
- nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.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<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);
-
- Dune::FieldVector<RF, dim> xgradient(0.0);
- xgradient[0] = 1000*0.353/0.684;
-
- RF rho_f=param.phys.rho_1*ch.pf;
- RF rho_f_reg= rho_f + param.delta;
-
- 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 mu_reg = pf_f*param.phys.mu;
-
- 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(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<ns.vfemspace.size(); i++)
- {
- for(int k=0; k<dim; k++)
- {
- //SIMPLIFY S=0, TODO: Implement this correctly
- 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 * ((ns.gradp[k] -xgradient[k])* ns.vbasis.phi[i]) //Pressure in strong form
- + 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
- + mu_reg * 8 / (param.CellHeight * param.CellHeight) * ns.v[k]*ns.vbasis.phi[i] //Drag from bottom and top plate
- )*factor);
- // integrate (grad u)^T * grad phi_i //Transpose for Viscosity
- for(int dd=0; dd<dim; ++dd)
- r.accumulate(child(ns.vspace,k),i, param.phys.mu * (ns.jacv[dd][k] * (ns.vbasis.gradphi[i][0][dd])) * factor);
- }
- }
- }
- }
-
- };
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
index 3c4a2dac985ee072b0d657b6b0cb1841fddea668..fcbea52d26293335a7c03cfea56f90652d51afa8 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
@@ -40,7 +40,6 @@ private:
public:
enum {doPatternVolume=true};
enum {doAlphaVolume=true};
- enum { doAlphaBoundary = true };
//constructor
FEM_FS_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
@@ -56,7 +55,7 @@ public:
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
+ const int order = 5;
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
@@ -82,16 +81,16 @@ public:
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 Rspeed = param.dw.eval_q(ch.pfOld);
RF R1 = -Rspeed/param.eps*
- (param.reaction.eval(ch.u)+param.phys.curvatureAlpha*ch.xi)*ch.basis.phi[i]; //TODO Scale with sigma13
+ (param.reaction.eval(ch.u)+param.phys.curvatureAlpha*ch.xi)*ch.basis.phi[i];
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(ch.pfspace,i,factor*(
(ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+J1
- //+(pf1*div_v + (gradpf1*v))*phi[i] TODO ADD THIS AGAIN WHEN BC ARE FIXED
+ +(pf1*div_v + (gradpf1*v))*phi[i]
-R1
));
@@ -118,57 +117,3 @@ public:
}
}
}
-
- // 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 e991d24fbba3d0c50f7e9a172cf98b66e0599001..1037e8ba16dedc3313cb2bbd57412dab198ce747 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
@@ -39,21 +39,12 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
public:
- //! Boundary condition indicator type
- // using BC = StokesBoundaryCondition; // ??
-
-
-// static const bool navier = P::assemble_navier;
-// static const bool full_tensor = P::assemble_full_tensor;
// pattern assembly flags
enum { doPatternVolume = true };
// residual assembly flags
enum { doAlphaVolume = true };
- //enum { doLambdaVolume = true };
-
- //TODO enum { doLambdaBoundary = true };
FEM_FS_CHNS_R_InstatStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
@@ -113,7 +104,6 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
+ 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)
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh b/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh
index 43d85cf860530a899871c16ac44b920cc782cba8..46fef46ac3a79db1b038df9b59d3a4bd4f75333e 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh
@@ -39,21 +39,13 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
public:
- //! Boundary condition indicator type
- // using BC = StokesBoundaryCondition; // ??
-// static const bool navier = P::assemble_navier;
-// static const bool full_tensor = P::assemble_full_tensor;
-
// pattern assembly flags
enum { doPatternVolume = true };
// residual assembly flags
enum { doAlphaVolume = true };
- //enum { doLambdaVolume = true };
-
- //TODO enum { doLambdaBoundary = true };
FEM_FS_CHNS_R_ThinPlateFlow (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
@@ -71,7 +63,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
+ const int order = 5;
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chLocal.DoEvaluation(eg);
@@ -105,15 +97,12 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
{
for(int k=0; k<dim; k++)
{
- //SIMPLIFY S=0, TODO: Implement this correctly
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 * (ns.gradp[k] * ns.vbasis.phi[i]) //Pressure in strong form
+ 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
+ mu_reg * 8 / (param.CellHeight * param.CellHeight) * ns.v[k]*ns.vbasis.phi[i] //Drag from bottom and top plate
)*factor);
// integrate (grad u)^T * grad phi_i //Transpose for Viscosity
diff --git a/dune/phasefield/porenetwork/1p_chns.hh b/dune/phasefield/porenetwork/1p_chns.hh
deleted file mode 100644
index d05c13d6099378ae5476c2766c6de37921357dc6..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/1p_chns.hh
+++ /dev/null
@@ -1,183 +0,0 @@
-#pragma once
-
-
-#include<dune/pdelab/constraints/common/constraints.hh>
-#include<dune/pdelab/constraints/conforming.hh>
-
-#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
-#include<dune/pdelab/gridfunctionspace/subspace.hh>
-#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
-
-//_________________________________________________________________________________________________________________________
-
-
-template<typename RF, typename ES, typename V_FEM, typename PHI_FEM>
-class GFS_1p_chns
-{
-public:
- typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
- typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
- typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
- typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
- typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
-
- ConstraintsAssembler conphi, conmu;
-
- typedef Dune::PDELab::GridFunctionSpace<ES, V_FEM, ConstraintsAssembler, VectorBackend> VX_GFS;
- VX_GFS vxGfs;
-
- typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> PHI_GFS;
- PHI_GFS phiGfs;
-
- typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> MU_GFS;
- MU_GFS muGfs;
-
- typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, VX_GFS> NS;
- NS ns;
-
- typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH;
- CH ch;
-
- typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
- NS_CC nsCC;
-
- typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
- CH_CC chCC;
-
- GFS_1p_chns(const ES& es, const V_FEM& vFem, const PHI_FEM& phiFem) :
- conphi(),conmu(),
- vxGfs(es,vFem),
- phiGfs(es,phiFem, conphi),
- muGfs(es,phiFem, conmu),
- ns(vxGfs),
- ch(phiGfs, muGfs),
- nsCC(),
- chCC()
- {
- vxGfs.name("velocity");
- phiGfs.name("phasefield");
- muGfs.name("potential");
- }
-
- template<typename Bdry>
- void updateConstraintsContainer(const Bdry& bdry)
- {
- Dune::PDELab::constraints(bdry.bdryNS,ns,nsCC, 0);
- Dune::PDELab::constraints(bdry.bdryCH,ch,chCC, 0);
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
-
-template<typename GFS, typename NS_V, typename CH_V>
-class DGF_1p_chns
-{
-public:
-
- typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VXSubGFS;
- VXSubGFS vxSubGfs;
- typedef Dune::PDELab::DiscreteGridFunction<VXSubGFS,NS_V> VXDGF;
- VXDGF vxDgf;
-
- typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<0> > PhiSubGFS;
- PhiSubGFS phiSubGfs;
- typedef Dune::PDELab::DiscreteGridFunction<PhiSubGFS,CH_V> PhiDGF;
- PhiDGF phiDgf;
-
- typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<1> > MuSubGFS;
- MuSubGFS muSubGfs;
- typedef Dune::PDELab::DiscreteGridFunction<MuSubGFS,CH_V> MuDGF;
- MuDGF muDgf;
-
-
- DGF_1p_chns(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
- vxSubGfs(gfs.ns),
- vxDgf(vxSubGfs,nsV),
- phiSubGfs(gfs.ch),
- phiDgf(phiSubGfs,chV),
- muSubGfs(gfs.ch),
- muDgf(muSubGfs,chV)
- {
- }
-
- template<typename VTKWRITER>
- void addToVTKWriter(VTKWRITER& vtkwriter)
- {
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VXDGF> >(vxDgf,"vx"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
-template<typename IniVX, typename IniPhi, typename IniMu>
-class Ini_1p_chns
-{
-public:
- IniVX iniVx;
- IniPhi iniPhi;
- IniMu iniMu;
- typedef Dune::PDELab::CompositeGridFunction<IniVX> IniNS;
- typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu> IniCH;
- IniNS iniNS;
- IniCH iniCH;
-
- template<typename GV, typename Params>
- Ini_1p_chns(const GV& gv, const Params& params) :
- iniVx(gv, params),
- iniPhi(gv, params),
- iniMu(gv, params),
- iniNS(iniVx),
- iniCH(iniPhi,iniMu)
- {
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
-template<typename BdryVX, typename BdryPhi, typename BdryMu>
-class Bdry_1p_chns
-{
-public:
- BdryVX bdryVx;
- BdryPhi bdryPhi;
- BdryMu bdryMu;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNS;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
- BdryNS bdryNS;
- BdryCH bdryCH;
-
- Bdry_1p_chns() :
- bdryVx(),
- bdryPhi(),
- bdryMu(),
- bdryNS(bdryVx),
- bdryCH(bdryPhi,bdryMu)
- {
- }
-};
-
-template<typename BdryVX, typename BdryPhi, typename BdryMu>
-class Bdry_1p_chns_Explicit_Domain
-{
-public:
- BdryVX bdryVx;
- BdryPhi bdryPhi;
- BdryMu bdryMu;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNS;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
- BdryNS bdryNS;
- BdryCH bdryCH;
-
- template<typename DomainSize>
- Bdry_1p_chns_Explicit_Domain(DomainSize ds) :
- bdryVx(ds),
- bdryPhi(ds),
- bdryMu(ds),
- bdryNS(bdryVx),
- bdryCH(bdryPhi,bdryMu)
- {
- }
-};
diff --git a/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh b/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
deleted file mode 100644
index b34ec3ef4a5e329798f6da4d5bb9c3ff57416be8..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
+++ /dev/null
@@ -1,88 +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>
-class FEM_1P_CahnHilliard_R :
- public Dune::PDELab::FullVolumePattern,
- public Dune::PDELab::LocalOperatorDefaultFlags,
- public Dune::PDELab::NumericalJacobianVolume<FEM_1P_CahnHilliard_R<Param, CHGFS, CHContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<FEM_1P_CahnHilliard_R<Param, CHGFS, CHContainer> >
-{
-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_1P_CahnHilliard_R(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) :
- 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_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 Rspeed = ch.pfOld < 0.95 ? param.dw.eval_q(ch.pfOld) : 0.0f;
- RF R1 = -Rspeed/param.eps*param.reaction.eval(1)*ch.basis.phi[i];
- //std::cout << "mobility" << param.mobility.eval(ch.pf) << "R1" << R1 << "dt" << param.time.dt << std::endl;
- RF J1 = param.mobility.eval(ch.pf) * (ch.gradxi*ch.basis.gradphi[i][0]);
-
- // 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
- +J1-R1
- ));
-
- 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])
- ));
- }
- }
- }
-
-};
diff --git a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
deleted file mode 100644
index aa3294a31c334a205aadb3375e4bbdf44043706f..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
+++ /dev/null
@@ -1,91 +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"fem_base_porenetwork.hh"
-
-
-
-
-template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
- class FEM_1P_NavierStokes :
- public Dune::PDELab::FullVolumePattern,
- public Dune::PDELab::LocalOperatorDefaultFlags,
- public Dune::PDELab::NumericalJacobianVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
- {
- private:
- typedef typename CHContainer::field_type RF;
- typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
-
- CHLocalBasisCache<CHGFS> chCache;
- Micro_VXLocalBasisCache<NSGFS> vCache;
-
- typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
- mutable CHLView chLocal;
-
- Param& param; // parameter functions
-
- //FOR DETECTING ENTITY changes
- typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
- mutable ENT OldEnt;
-
- public:
-
-
- // pattern assembly flags
- enum { doPatternVolume = true };
-
- // residual assembly flags
- enum { doAlphaVolume = true };
- //enum { doLambdaVolume = true };
-
-
-
- FEM_1P_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
- NSGFS& nsGfs_) :
- chLocal(chGfs_,chV_),
- 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
- {
-// std::cout << "hi_alpha" << std::endl;
- //Quadrature Rule
- const int dim = EG::Entity::dimension;
- const int order = 5; //TODO: Choose more educated
- auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
-
- chLocal.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<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
- NSPorenetworkVariables<RF,dim,LFSV> ns(vCache, ip.position(), lfsv, x.base(), S);
- RF pf_f=ch.pf+param.delta;
-
- for (size_t i=0; i<ns.vxspace.size(); i++) //TODO Check
- {
- r.accumulate(ns.vxspace, i, factor*(
- param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]
- + param.phys.eval_mu(pf_f) * (ns.gradvx * ns.vbasis.gradphi[i][0])
- - pf_f * ns.vbasis.phi[i]
- ));
- } // for i
- } // for ip
-// std::cout << "bye_alpha" << std::endl;
- } //alpha volume
-
- };
diff --git a/dune/phasefield/porenetwork/fem_base_porenetwork.hh b/dune/phasefield/porenetwork/fem_base_porenetwork.hh
deleted file mode 100644
index ceb0b62d253c5ef93ce2b22dd24630039f8cf301..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/fem_base_porenetwork.hh
+++ /dev/null
@@ -1,34 +0,0 @@
-#pragma once
-
-#include<dune/phasefield/localoperator/fem_base_chns.hh>
-
-
-template<typename RF, const int dim, typename NSLFS>
-class NSPorenetworkVariables
-{
-public:
- const typename NSLFS::template Child<0>::Type& vxspace;
- BasisFN<typename NSLFS::template Child<0>::Type> vbasis;
-
- RF vx;
- Dune::FieldVector<RF,dim> gradvx;
-
- template<typename Pos, typename vCache, typename LV, typename S>
- NSPorenetworkVariables(const vCache& vcache, const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
- vxspace(child(nsLfs,Dune::Indices::_0)),
- vbasis(vcache, position, vxspace, s),
- vx(0.0),
- gradvx(0.0)
- {
- for (size_t i=0; i<this->vxspace.size(); i++)
- {
- vx += lv[vxspace.localIndex(i)]*vbasis.phi[i];
- gradvx.axpy(lv[vxspace.localIndex(i)],vbasis.gradphi[i][0]);
- }
- }
-};
-
-
-template<typename NSGFS>
-using Micro_VXLocalBasisCache =
- Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
deleted file mode 100644
index fa1e81f815e7df82dd2a134374ecad2b97642e27..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/pn_1pflow.hh
+++ /dev/null
@@ -1,292 +0,0 @@
-#pragma once
-
-// C includes
-#include<sys/stat.h>
-// C++ includes
-#include<iostream>
-#include<vector>
-// Dune includes
-#include<dune/grid/io/file/vtk.hh>
-#include<dune/pdelab/newton/newton.hh>
-#include<dune/pdelab/common/functionutilities.hh>
-#include<dune/pdelab/function/callableadapter.hh>
-
-#include <dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh>
-#include <dune/phasefield/porenetwork/fem_1p_navierstokes.hh>
-#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
-
-#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
-
-#include <dune/phasefield/boundary_rectangular.hh>
-#include <dune/phasefield/porenetwork/pn_initial.hh>
-#include <dune/phasefield/chns_base.hh>
-#include <dune/phasefield/porenetwork/1p_chns.hh>
-#include <dune/phasefield/initial_fn.hh>
-#include <dune/phasefield/debug/vector_debug.hh>
-
-template< typename Grid>
-std::shared_ptr<Grid> initPnGrid(std::string filename)
-{
- Dune::GridFactory<Grid> factory;
- Dune::GmshReader<Grid>::read(factory,filename,true,true);
- return std::shared_ptr<Grid>(factory.createGrid());
-}
-
-
-template<typename Parameters>
-class Pn_1PFlow{
-private:
-
- typedef double RF;
- const int dim = 2;
-
- typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
- typedef Grid::ctype DF;
- typedef Grid::LeafGridView GV;
- typedef Dune::PDELab::AllEntitySet<GV> ES;
- typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,2> V_FEM;
- typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> PHI_FEM;
-
- typedef GFS_1p_chns<RF, ES, V_FEM, PHI_FEM> GFS;
-
- typedef typename GFS::NS NS_GFS;
- typedef typename GFS::CH CH_GFS;
- typedef Dune::PDELab::Backend::Vector<CH_GFS, RF> CH_V;
- typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
-
- typedef Ini_1p_chns< ZeroInitial<GV,RF,Parameters>, //V
- PhiInitial_Quartersquare<GV,RF,Parameters>,
- ZeroInitial<GV,RF,Parameters>> Initial;
- typedef DomainSize<RF> DS;
- typedef Bdry_1p_chns_Explicit_Domain<LeftBottomBoundary<DS>, //V
- NoBoundary<DS>,
- NoBoundary<DS>> Bdry;
-
- typedef FEM_1P_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
- typedef FEM_1P_CahnHilliard_R<Parameters, CH_GFS, CH_V> CH_LOP;
- typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
-
- typedef DGF_1p_chns<GFS, NS_V, CH_V> DGF;
-
- typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
- typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
-
- typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS;
- typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS;
-
- typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> NS_Newton;
- typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> CH_Newton;
-
- Parameters& param;
- std::shared_ptr<Grid> gridp;
- GV gv;
- ES es;
- V_FEM vFem;
- PHI_FEM phiFem;
-
- GFS gfs;
-
- CH_V chV, chVOld;
- NS_V nsV;
- Initial initial;
- DS ds;
- Bdry bdry;
-
- MBE mbe;
-
- RF kf = 0;
- RF phiSolid = 0;
-
- std::unique_ptr<NS_LOP> nsLop;
- std::unique_ptr<CH_LOP> chLop;
- std::unique_ptr<DGF> dgf;
- std::unique_ptr<NS_GO> nsGo;
- std::unique_ptr<CH_GO> chGo;
- std::unique_ptr<NS_LS> nsLs;
- std::unique_ptr<CH_LS> chLs;
- std::unique_ptr<NS_Newton> nsNewton;
- std::unique_ptr<CH_Newton> chNewton;
- std::unique_ptr<Dune::VTKSequenceWriter<GV>> vtkwriter;
-
- //------------------------------------------------------------------------------------
-
-
- void initializeGridOperators()
- {
- //Grid Operator
- nsGo = std::make_unique<NS_GO>(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nsLop,mbe);
- chGo = std::make_unique<CH_GO>(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,*chLop,mbe);
- // Linear solver
- chLs = std::make_unique<CH_LS>(*chGo,500,100,3,1);
- nsLs = std::make_unique<NS_LS>(*nsGo,1000,200,3,1);
- // Nonlinear solver
- nsNewton = std::make_unique<NS_Newton>(*nsGo,nsV,*nsLs);
- nsNewton->setParameters(param.nsNewtonTree);
-
- chNewton = std::make_unique<CH_Newton>(*chGo,chV,*chLs);
- chNewton->setParameters(param.chNewtonTree);
-
- //printVector(chV,10,"chV");
- //printVector(nsV,10,"nsV");
- }
-
- void remesh()
- {
- // mark elements for refinement
- typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
- RIND rInd(gv,gfs.ch,param.pTree.sub("Refinement"));
- rInd.calculateIndicator(gv,chV,false);
- rInd.markGrid(*gridp,2);
- rInd.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
- // recompute constraints
- gfs.updateConstraintsContainer(bdry);
- interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
- initializeGridOperators();
- }
-
- void applyNs(int verbosity = 0)
- {
- if (verbosity > 0) std::cout << "= Begin newtonapply NS:" << std::endl;
- nsNewton->apply();
- }
-
- void applyCh(int verbosity = 0)
- {
- chVOld = chV;
- if (verbosity > 0) std::cout << "= Begin newtonapply CH:" << std::endl;
- chNewton->apply();
- //chVOld = chV;
- remesh();
- }
-
- void updateKf()
- {
- auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
- Dune::FieldVector<double,1> w;
- dgf->vxDgf.evaluate(element, xlocal, w);
- Dune::FieldVector<double,1> phi1;
- dgf->phiDgf.evaluate(element, xlocal, phi1);
- return Dune::FieldVector<double,1>(phi1*w);
- //Dune::FieldVector<double,1> phi2;
- //dgf->phi2Dgf.evaluate(element, xlocal, phi2);
- //return Dune::FieldVector<double,1>((phi1+phi2)*w);
- });
- Dune::FieldVector<double,1> integral;
- Dune::PDELab::integrateGridFunction(productFunction,integral,2);
- kf = integral[0];
- }
-
- void updatePhiSolid()
- {
- auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
- Dune::FieldVector<double,1> phi1;
- dgf->phiDgf.evaluate(element, xlocal, phi1);
- return Dune::FieldVector<double,1>(1-phi1);
- //Dune::FieldVector<double,1> phi2;
- //dgf->phi2Dgf.evaluate(element, xlocal, phi2);
- //return Dune::FieldVector<double,1>(1-phi1-phi2);
- });
- Dune::FieldVector<double,1> integral;
- Dune::PDELab::integrateGridFunction(productFunction,integral,2);
- phiSolid = integral[0];
- }
-
-//------------------------------------------------------------------------------------
-
-public:
-
- //Constructor
- Pn_1PFlow(std::string GridFilename, std::string OutputFilename, Parameters& param_) :
- param(param_),
- gridp(initPnGrid<Grid>(GridFilename)),
- gv(gridp->leafGridView()),
- es(gv),
- vFem(es), phiFem(es),
- gfs(es,vFem,phiFem),
- chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch),
- initial(gv,param), ds(0,1,0,1,1e-6), bdry(ds),mbe(10)
- {
- Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
- Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
-
- // Refine Intial Values
- for(int i=0; i<param.pTree.get("Refinement.maxRefineLevel",2); i++)
- {
- // mark elements for refinement
- typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND0;
- RIND0 rInd0(gv,gfs.ch,param.pTree.sub("Refinement"));
- rInd0.calculateIndicator(gv,chV,false);
- rInd0.markGrid(*gridp,2);
- rInd0.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
- //rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);
- Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
- Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
- }
-
- gfs.updateConstraintsContainer(bdry);
- interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
-
- //chVOld = chV;
-
- //Local Operator
- nsLop = std::make_unique<NS_LOP>(param, gfs.ch, chV, gfs.ns);
- chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld);
-
- dgf = std::make_unique<DGF>(gfs,nsV,chV);
-
- struct stat st;
- if( stat( OutputFilename.c_str(), &st ) != 0 )
- {
- int stat = 0;
- stat = mkdir(OutputFilename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
- if( stat != 0 && stat != -1)
- std::cout << "Error: Cannot create directory " << OutputFilename << std::endl;
- }
- auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
- vtkwriter = std::make_unique<Dune::VTKSequenceWriter<GV> >(stationaryvtkwriter,OutputFilename.c_str(),OutputFilename.c_str(),"");
- dgf->addToVTKWriter(*vtkwriter);
-
-
- initializeGridOperators();
- applyNs();
- updatePhiSolid();
- updateKf();
- }
-
- void writeVTK(RF time)
- {
- vtkwriter->write(time,Dune::VTK::ascii);
- }
-
- void grow(RF length)
- {
- if(length <= 0)
- {
- std::cout << "Need positive length" << std::endl;
- return;
- }
- RF dtmax = param.pTree.get("Time.dtMax",0.01);
- while(length > dtmax)
- {
- writeVTK(1-length);
- length -= dtmax;
- param.time.dt = dtmax;
- applyCh(0);
- }
- param.time.dt = length;
- applyCh(0);
- applyNs(0);
- updateKf();
- updatePhiSolid();
- }
-
- RF getKf()
- {
- return kf;
- }
-
- RF getPhiSolid()
- {
- return phiSolid;
- }
-};
diff --git a/dune/phasefield/porenetwork/pn_initial.hh b/dune/phasefield/porenetwork/pn_initial.hh
deleted file mode 100644
index df0c6218efe1f24cbcbf031e0b6809713968ee8e..0000000000000000000000000000000000000000
--- a/dune/phasefield/porenetwork/pn_initial.hh
+++ /dev/null
@@ -1,33 +0,0 @@
-#pragma once
-
-template<typename GV, typename RF, typename Params>
-class PhiInitial_Quartersquare :
- public Dune::PDELab::AnalyticGridFunctionBase<
- Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Quartersquare<GV,RF,Params> >
-{
- const Params& param;
- RF eps;
- RF r;
-
-public:
- typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
-
- PhiInitial_Quartersquare(const GV & gv, const Params & params) :
- Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Quartersquare<GV,RF,Params> > (gv),
- param(params)
- {
- eps = params.eps;
- r = params.initial.ptree.get("Radius",0.2);
- std::cout << "EPS" << eps;
- }
-
- inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
- {
- RF rf1 = param.dw.eval_shape((x[0]-r)*1/eps);
- RF rf2 = param.dw.eval_shape((x[1]-r)*1/eps);
- //RF rf1 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-r)));
- //RF rf2 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-r)));
-
- y=rf1*rf2;
- }
-};