From 2a58ccabc0c5b4bd40981ef577e12d23ec0e1266 Mon Sep 17 00:00:00 2001
From: Lars von Wolff <lars.von-wolff@ians.uni-stuttgart.de>
Date: Fri, 30 Aug 2019 15:12:00 +0200
Subject: [PATCH] +added missing files of previous commit
---
dune/phasefield/{ => archive}/acfs_ibvp.hh | 0
dune/phasefield/{ => archive}/acfs_problem.hh | 0
dune/phasefield/archive/ff_parameters.hh | 160 ++++++
.../tmp___myincludes.hh} | 0
dune/phasefield/base_parameters.hh | 515 ++++++++++++++++++
dune/phasefield/boundary_fn.hh | 71 +++
dune/phasefield/chns_base.hh | 263 ++-------
dune/phasefield/debug/vector_debug.hh | 2 +-
dune/phasefield/ff_chns.hh | 181 ++++++
dune/phasefield/ffs_chns_r.hh | 232 ++++++++
dune/phasefield/fs_chns_r.hh | 198 +++++++
dune/phasefield/initial_fn.hh | 162 +++++-
.../phasefield/localoperator/fem_base_chns.hh | 298 ++++++++++
...illiard.hh => fem_ff_chns_cahnhilliard.hh} | 22 +-
...rstokes.hh => fem_ff_chns_navierstokes.hh} | 28 +-
.../fem_ffs_chns_navierstokes.hh | 239 ++++++++
.../fem_ffs_chns_r_cahnhilliard.hh | 226 ++++++++
.../fem_fs_chns_r_cahnhilliard.hh | 217 ++++++++
.../fem_fs_chns_r_instatstokes.hh | 246 +++++++++
.../localoperator/pf_3p_diff_estimator.hh | 84 +++
.../localoperator/pf_diff_estimator.hh | 13 +-
.../ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh | 12 +-
dune/phasefield/utility/pf_3p_color_dgf.hh | 112 ++++
src/CMakeLists.txt | 2 +
src/test_2f0s/{chns.ini => 2f0s.ini} | 49 +-
src/test_2f0s/CMakeLists.txt | 2 +-
src/test_2f0s/chns_driver.hh | 122 ++---
src/test_2f0s/chns_equation.cc | 4 +-
src/test_3pchns_r/3p_driver.hh | 216 ++++++++
src/test_3pchns_r/3p_equation.cc | 116 ++++
src/test_3pchns_r/CMakeLists.txt | 2 +
src/test_3pchns_r/chns.ini | 58 ++
src/test_3pchns_r/myincludes.hh | 59 ++
33 files changed, 3540 insertions(+), 371 deletions(-)
rename dune/phasefield/{ => archive}/acfs_ibvp.hh (100%)
rename dune/phasefield/{ => archive}/acfs_problem.hh (100%)
create mode 100644 dune/phasefield/archive/ff_parameters.hh
rename dune/phasefield/{myincludes.hh => archive/tmp___myincludes.hh} (100%)
create mode 100644 dune/phasefield/base_parameters.hh
create mode 100644 dune/phasefield/boundary_fn.hh
create mode 100644 dune/phasefield/ff_chns.hh
create mode 100644 dune/phasefield/ffs_chns_r.hh
create mode 100644 dune/phasefield/fs_chns_r.hh
create mode 100644 dune/phasefield/localoperator/fem_base_chns.hh
rename dune/phasefield/localoperator/{chns_fem_cahnhilliard.hh => fem_ff_chns_cahnhilliard.hh} (85%)
rename dune/phasefield/localoperator/{chns_fem_navierstokes.hh => fem_ff_chns_navierstokes.hh} (87%)
create mode 100644 dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
create mode 100644 dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
create mode 100644 dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
create mode 100644 dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
create mode 100644 dune/phasefield/localoperator/pf_3p_diff_estimator.hh
rename src/test_2f0s/acfs_fem_ac_estimator.hh => dune/phasefield/localoperator/pf_diff_estimator.hh (82%)
create mode 100644 dune/phasefield/utility/pf_3p_color_dgf.hh
rename src/test_2f0s/{chns.ini => 2f0s.ini} (59%)
create mode 100644 src/test_3pchns_r/3p_driver.hh
create mode 100644 src/test_3pchns_r/3p_equation.cc
create mode 100644 src/test_3pchns_r/CMakeLists.txt
create mode 100644 src/test_3pchns_r/chns.ini
create mode 100644 src/test_3pchns_r/myincludes.hh
diff --git a/dune/phasefield/acfs_ibvp.hh b/dune/phasefield/archive/acfs_ibvp.hh
similarity index 100%
rename from dune/phasefield/acfs_ibvp.hh
rename to dune/phasefield/archive/acfs_ibvp.hh
diff --git a/dune/phasefield/acfs_problem.hh b/dune/phasefield/archive/acfs_problem.hh
similarity index 100%
rename from dune/phasefield/acfs_problem.hh
rename to dune/phasefield/archive/acfs_problem.hh
diff --git a/dune/phasefield/archive/ff_parameters.hh b/dune/phasefield/archive/ff_parameters.hh
new file mode 100644
index 0000000..1d6a836
--- /dev/null
+++ b/dune/phasefield/archive/ff_parameters.hh
@@ -0,0 +1,160 @@
+
+#pragma once
+//#include "../myincludes.hh"
+#include <cmath>
+
+template<typename Number>
+class Problem
+{
+ Number t;
+
+ const Number TurningPointLeft = 0.5 - 1/sqrt(12);
+ const Number TurningPointRight = 0.5 + 1/sqrt(12);
+public:
+ typedef Number value_type;
+ Number eps, delta;
+ Number D;
+ Number dt;
+ Number uAst;
+ Number u_regularization;
+ Number mu;
+
+ Number rho_0;
+
+ Number rho_1;
+ Number rho_2;
+
+ Number M;
+ Number M_reduced;
+ Number fscaling;
+ Number SurfaceTensionScaling;
+
+ //Number rho (const Number& x) {return rho_0;}
+ Number rho (const Number& x){return x*rho_2 + (1-x)*rho_1;}
+
+ Number Sigma1,Sigma2,Sigma3,SigmaT;
+
+
+ // Constructor takes parameter
+ Problem (const Number& dt_, const Dune::ParameterTree& ptree)
+ : t(0.0)
+ {
+ eps = ptree.get("problem.eps",(Number)0.1);
+ delta = ptree.get("problem.delta",(Number)0.1);
+ M = ptree.get("problem.Mobility",(Number)0.1);
+ M_reduced = ptree.get("problem.MobilityReduced",(Number)0.1);
+ u_regularization = ptree.get("problem.URegularization",(Number)0.01);
+
+ D = ptree.get("physics.D",(Number)0.1);
+ uAst = ptree.get("physics.uAst",(Number)0.1);
+ rho_0 = ptree.get("physics.rho",(Number)1);
+ mu = ptree.get("physics.mu",(Number)0.001);
+ fscaling = ptree.get("phsyics.fscaling",(Number)1);
+ SurfaceTensionScaling = ptree.get("physics.SurfaceTensionScaling",(Number)1);
+ rho_1 = ptree.get("3p.rho1",(Number)1);
+ rho_2 = ptree.get("3p.rho2",(Number)1);
+
+ Sigma1 = ptree.get("3p.Sigma1",(Number)1);
+ Sigma2 = ptree.get("3p.Sigma2",(Number)1);
+ Sigma3 = ptree.get("3p.Sigma3",(Number)1);
+ SigmaT = 1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3);
+
+ dt=dt_;
+ }
+
+ Number Mob (const Number& u) const
+ {
+ return M; //M_reduced + (M-M_reduced)*q(u);
+ }
+
+ //Double Well
+ Number W (const Number& u) const
+ {
+ return u*u*(1-u)*(1-u);
+ }
+ // Double Well Derivative(/Secant)
+ Number dW (const Number& x, const Number& y) const
+ {
+ return (x+y)*(1 +x*x +y*y) -2*(x*x +x*y +y*y);
+ }
+
+ // Double Well Derivative
+ Number dW (const Number& x) const
+ {
+ return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
+ }
+
+ // Double Well Derivative of Convex Part
+ Number dWConvex (const Number& x) const
+ {
+ if(x<TurningPointLeft)
+ return dW(x)-dW(TurningPointLeft);
+ if(x>TurningPointRight)
+ return dW(x)-dW(TurningPointRight);
+ return 0;
+ }
+
+ // Double Well Derivative of Convcarve Part
+ Number dWConcarve (const Number& x) const
+ {
+ if(x<TurningPointLeft)
+ return dW(TurningPointLeft);
+ if(x>TurningPointRight)
+ return dW(TurningPointRight);
+ return dW(x);
+ }
+
+ Number dWdpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma1-SigmaT)*dW(pf1)-SigmaT*dW(pf2)-SigmaT*dW(pf3);
+ }
+
+ Number dWdpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma2-SigmaT)*dW(pf2)-SigmaT*dW(pf1)-SigmaT*dW(pf3);
+ }
+
+ Number dWdpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (Sigma3-SigmaT)*dW(pf3)-SigmaT*dW(pf1)-SigmaT*dW(pf2);
+ }
+
+ Number f (const Number& u) const
+ {
+ if(u<0)
+ return fscaling*(-1.0+1000*u);
+ else
+ return fscaling*(u*u-1);
+ }
+
+ Number q (const Number& pf) const
+ {
+ if (pf<0 || pf>1)
+ return 0;
+ return pf*(1-pf);
+ }
+
+ Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ if (pf1<0 || pf3<0)
+ return 0;
+ return pf1*pf3;
+ }
+
+ Number Diffusion(const Number& pf_c)
+ {
+ return D*fmax(1e-3,(pf_c-0.05)/0.95);
+ }
+
+ Number VelDissipation(const Number& pf_f)
+ {
+ return 10*(1-pf_f)*(1-pf_f);
+ }
+
+ //! Set time in instationary case
+ void setTime (Number t_)
+ {
+ t = t_;
+ }
+
+};
diff --git a/dune/phasefield/myincludes.hh b/dune/phasefield/archive/tmp___myincludes.hh
similarity index 100%
rename from dune/phasefield/myincludes.hh
rename to dune/phasefield/archive/tmp___myincludes.hh
diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
new file mode 100644
index 0000000..2089d4e
--- /dev/null
+++ b/dune/phasefield/base_parameters.hh
@@ -0,0 +1,515 @@
+
+#pragma once
+//#include "../myincludes.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>
+class Mobility_Const
+{
+public:
+ const Number mobility;
+ Mobility_Const(const Dune::ParameterTree & param) :
+ mobility(param.get("Mobility",(Number)0.1))
+ {
+ }
+
+ Number eval( const Number& pf) const
+ {
+ return mobility;
+ }
+
+ Number eval( const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return mobility;
+ }
+};
+
+template<typename Number>
+class Mobility_Degen
+{
+public:
+ const Number mobility;
+ const Number mobility_reg;
+ Mobility_Degen(const Dune::ParameterTree & param) :
+ mobility(param.get("Mobility",(Number)0.1)),
+ mobility_reg(param.get("MobilityMin",(Number)0.0))
+ {
+ }
+
+ Number eval( const Number& pf) const
+ {
+ return std::max(mobility*pf*(1-pf),mobility_reg);
+ }
+
+ Number eval( const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ std::cout << "Mobility_Degen for three phases not properly implemented yet" << std::endl;
+ return mobility;
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template<typename Number>
+class Reaction_Quadratic_Limited
+{
+public:
+ const Number reactionRate;
+ const Number reactionLimiter;
+ Reaction_Quadratic_Limited(const Dune::ParameterTree & param) :
+ reactionRate(param.get("ReactionRate",(Number)0.1)),
+ reactionLimiter(param.get("ReactionLimiter",(Number)100))
+ {
+ }
+
+ Number eval( const Number& u) const
+ {
+ if(u<0)
+ return reactionRate*(-1.0+reactionLimiter*u);
+ else
+ return reactionRate*(u*u-1.0);
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template<typename Number>
+class VelDissipation_Quadratic
+{
+public:
+ const Number velDissipationRate;
+ VelDissipation_Quadratic(const Dune::ParameterTree & param) :
+ velDissipationRate(param.get("VelDissipationRate",(Number)10))
+ {
+ }
+
+ Number eval( const Number& pf_f) const
+ {
+ return velDissipationRate*(1-pf_f)*(1-pf_f);
+ }
+};
+
+//______________________________________________________________________________________________________
+
+
+template<typename Number>
+class Params_Physical
+{
+public:
+ const Number D;
+ const Number uAst;
+ const Number mu;
+ const Number rho_1;
+ const Number rho_2;
+ const Number rho_3;
+ const Number curvatureAlpha;
+
+ Params_Physical(const Dune::ParameterTree& ptree, const Number eps) :
+ D(ptree.get("D",(Number)0.1)),
+ uAst(ptree.get("uAst",(Number)2)),
+ mu(ptree.get("mu",(Number)0.01)),
+ rho_1(ptree.get("rho1",(Number)1)),
+ rho_2(ptree.get("rho2",(Number)1)),
+ rho_3(ptree.get("rho1",(Number)1)), //WE ASSUME RHO1 = RHO3 !!!
+ curvatureAlpha(ptree.get("CurvatureAlpha", eps))
+ {
+ }
+
+ Number rho_f (const Number& pf1, const Number& pf2) const
+ {
+ return pf1*rho_1+pf2*rho_2;
+ }
+
+};
+
+//______________________________________________________________________________________________________
+class Params_Time
+{
+public:
+ double t;
+ double dt;
+ Params_Time(double t_, double dt_) :
+ t(t_),dt(dt_)
+ {
+ }
+};
+
+//______________________________________________________________________________________________________
+
+class Params_Initial
+{
+public:
+ const Dune::ParameterTree& ptree;
+ Params_Initial(const Dune::ParameterTree& tree) :
+ ptree(tree)
+ {
+ }
+};
+
+//______________________________________________________________________________________________________
+
+template <typename Number, typename Mobility>
+class Params_base
+{
+public:
+ Number eps;
+
+ Mobility mobility;
+ Params_Physical<Number> phys;
+ Params_Time time;
+
+ Params_Initial initial;
+
+
+ Params_base(const Dune::ParameterTree& ptree) :
+ eps(ptree.get("Phasefield.eps",(Number)0.1)),
+ mobility(ptree.sub("Phasefield")),
+ phys(ptree.sub("Physics"),eps),
+ time(0,0),
+ initial(ptree.sub("Initial"))
+ {
+ }
+};
+
+
+template <typename Number, typename DW, typename Mobility>
+class Params_ff : public Params_base<Number, Mobility>
+{
+public:
+
+ DW dw;
+ Number SurfaceTensionScaling;
+
+ Params_ff(const Dune::ParameterTree& ptree) :
+ Params_base<Number, Mobility>(ptree),
+ dw(ptree.sub("Phasefield")),
+ SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1))
+ {
+ }
+};
+
+
+template <typename Number, typename TW, typename Mobility, typename Reaction, typename VelDissipation>
+class Params_ffs_r : public Params_base<Number, Mobility>
+{
+public:
+ Number delta;
+
+ TW tw;
+ Reaction reaction;
+ VelDissipation vDis;
+
+ Number SurfaceTensionScaling;
+
+
+ Params_ffs_r(const Dune::ParameterTree& ptree) :
+ Params_base<Number, Mobility>(ptree),
+ delta(ptree.get("Phasefield.delta",this->eps)),
+ tw(ptree.sub("Phasefield")),
+ reaction(ptree.sub("Physics")),
+ vDis(ptree.sub("Physics")),
+ SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1))
+ {
+ }
+};
+
+
+template <typename Number, typename DW, typename Mobility, typename Reaction, typename VelDissipation>
+class Params_fs_r : public Params_base<Number, Mobility>
+{
+public:
+ Number delta;
+
+ DW dw;
+ Reaction reaction;
+ VelDissipation vDis;
+
+ Params_fs_r(const Dune::ParameterTree& ptree) :
+ Params_base<Number, Mobility>(ptree),
+ delta(ptree.get("Phasefield.delta",this->eps)),
+ dw(ptree.sub("Phasefield")),
+ reaction(ptree.sub("Physics")),
+ vDis(ptree.sub("Physics"))
+ {
+ }
+};
+//
+// 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/boundary_fn.hh b/dune/phasefield/boundary_fn.hh
new file mode 100644
index 0000000..06ea101
--- /dev/null
+++ b/dune/phasefield/boundary_fn.hh
@@ -0,0 +1,71 @@
+
+#pragma once
+#include <cmath>
+#include <dune/pdelab/common/function.hh>
+#include <dune/pdelab/constraints/common/constraintsparameters.hh>
+
+#define GRID_EPS 1e-7
+
+class InflowOutflow_v : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (xg[0]> 1-GRID_EPS && xg[1]>-1+GRID_EPS && xg[1]<1-GRID_EPS)
+ {
+ return false;
+ }
+ return true;
+ }
+};
+
+class InflowOutflow_p : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (xg[0]> 1-GRID_EPS && xg[1]>-1+GRID_EPS && xg[1]<1-GRID_EPS)
+ {
+ return true;
+ }
+ return false;
+ }
+};
+
+class InflowOutflow_u : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (xg[0]< GRID_EPS)
+ {
+ return true;
+ }
+ return false;
+ }
+};
+
+class InflowOutflow_phi : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (xg[0]< GRID_EPS || xg[0] >1-GRID_EPS)
+ {
+ return true;
+ }
+ return false;
+ }
+};
diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index 104a891..2405e7d 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -1,31 +1,13 @@
#pragma once
-#include "myincludes.hh"
-//#include <dune/xt/grid/grids.hh>
-//#include <dune/xt/grid/view/periodic.hh>
-#include <dune/pdelab/localoperator/stokesparameter.hh>
-#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
-#include <dune/pdelab/localoperator/taylorhoodnavierstokes.hh>
-#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
-#include <dune/pdelab/adaptivity/adaptivity.hh>
-
-//#include "acfs_fem_allencahn.hh"
-//#include "acfs_fem_cahnhilliard.hh"
-//#include "acfs_fem_navierstokes.hh"
-//#include "acfs_fem_ac_estimator.hh"
-//#include "threephase_fem_cahnhilliard.hh"
-//#include "threephase_fem_navierstokes.hh"
-//#include "threephase_dgf_coloring.hh"
-#include <dune/phasefield/localoperator/chns_fem_cahnhilliard.hh>
-#include <dune/phasefield/localoperator/chns_fem_navierstokes.hh>
-#include <dune/phasefield/acfs_problem.hh>
-#include <dune/phasefield/acfs_ibvp.hh>
-#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/pdelab/finiteelementmap/p0fem.hh>
+#include <dune/pdelab/adaptivity/adaptivity.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/backend/istl.hh>
-//template <class A>
-//using MySolver = Dune::RestartedGMResSolver<A, A, A>;
-
+#include <dune/phasefield/base_parameters.hh>
//_________________________________________________________________________________________________________________________
@@ -48,8 +30,8 @@ class TimestepManager
{
public:
double T;
- double dt;
- double t;
+
+ Params_Time& time;
double dtMin;
double dtMax;
@@ -61,21 +43,22 @@ public:
int verb;
- TimestepManager(const Dune::ParameterTree & param, int verbosity)
+ TimestepManager(const Dune::ParameterTree & param, Params_Time& time_, int verbosity) :
+ time(time_)
{
verb = verbosity;
T=param.get("tMax",(double)1.0);
- dt=param.get("dt",(double)0.1);
+ time.dt=param.get("dt",(double)0.1);
- dtMin = param.get("dtMin",(double)dt);
- dtMax = param.get("dtMax",(double)dt);
+ dtMin = param.get("dtMin",(double)time.dt);
+ dtMax = param.get("dtMax",(double)time.dt);
dtIncreaseFactor = param.get("dtIncreaseFactor",(double)1.2);
dtDecreaseFactor = param.get("dtDecreaseFactor",(double)0.5);
- saveintervall = param.get("temporal.saveintervall",(double)0.0);
+ saveintervall = param.get("saveintervall",(double)0.0);
- t=param.get("tIni",(double)0.0);
- lastsave = t - saveintervall;
+ time.t=param.get("tIni",(double)0.0);
+ lastsave = time.t - saveintervall;
}
// RETURN VALUE: -1 -> ABORT, 0-> REDO TIMESTEP, 1-> OK
@@ -84,15 +67,16 @@ public:
{
try
{
+ 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();
- if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
- ns_newton.apply();
}
catch(...)
{
- if(dt <= dtMin)
+ if(time.dt <= dtMin)
{
std::cout << "===============================" << std::endl;
if(dtMin==dtMax)
@@ -102,23 +86,23 @@ public:
std::cout << "===============================" << std::endl;
return -1;
}
- dt = std::max(dt*dtDecreaseFactor,dtMin);
- if (verb > 0) std::cout << "Reducing Timestep to "<< dt << std::endl;
+ time.dt = std::max(time.dt*dtDecreaseFactor,dtMin);
+ if (verb > 0) std::cout << "Reducing Timestep to "<< time.dt << std::endl;
return 0;
}
- t += dt;
- if(dt < dtMax)
+ time.t += time.dt;
+ if(time.dt < dtMax)
{
- dt = std::min(dt*dtIncreaseFactor, dtMax);
- if (verb > 0) std::cout << "New Timestep "<< dt << std::endl;
+ time.dt = std::min(time.dt*dtIncreaseFactor, dtMax);
+ if (verb > 0) std::cout << "New Timestep "<< time.dt << std::endl;
}
- if( t>lastsave+saveintervall )
+ if( time.t>lastsave+saveintervall )
{
- lastsave = t;
- vtkwriter.write(t,Dune::VTK::ascii);
+ lastsave = time.t;
+ vtkwriter.write(time.t,Dune::VTK::ascii);
}
return 1;
}
@@ -126,7 +110,7 @@ public:
bool isFinished()
{
- return (t>=T);
+ return (time.t>=T);
}
};
@@ -143,177 +127,6 @@ void interpolateToBdry(const GFS& gfs, const NS_BDRY_FN& nsBdryFn, const CH_BDRY
//_________________________________________________________________________________________________________________________
-template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
-class GFS_2f0s
-{
-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 conp, conphi, conmu;
-
- typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
- VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
- V_GFS vGfs;
-
- typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
- P_GFS pGfs;
-
- 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, V_GFS, P_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_2f0s(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
- conp(),conphi(),conmu(),
- vGfs(es,vFem),
- pGfs(es,pFem, conp),
- phiGfs(es,phiFem, conphi),
- muGfs(es,phiFem, conmu),
- ns(vGfs, pGfs),
- ch(phiGfs, muGfs),
- nsCC(),
- chCC()
- {
- vGfs.name("velocity");
- pGfs.name("pressure");
- 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_2f0s
-{
-public:
-
- typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
- VSubGFS vSubGfs;
- typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
- VDGF vDgf;
-
- typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
- PSubGFS pSubGfs;
- typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
- PDGF pDgf;
-
- 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_2f0s(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
- vSubGfs(gfs.ns),
- vDgf(vSubGfs,nsV),
- pSubGfs(gfs.ns),
- pDgf(pSubGfs,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<VDGF> >(vDgf,"v"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
-template<typename IniV, typename IniP, typename IniPhi, typename IniMu>
-class Ini_2f0s
-{
-public:
- IniV iniV;
- IniP iniP;
- IniPhi iniPhi;
- IniMu iniMu;
- typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
- typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu> IniCH;
- IniNS iniNS;
- IniCH iniCH;
-
- template<typename GV, typename Params>
- Ini_2f0s(const GV& gv, const Params& params) :
- iniV(gv, params),
- iniP(gv, params),
- iniPhi(gv, params),
- iniMu(gv, params),
- iniNS(iniV,iniP),
- iniCH(iniPhi,iniMu)
- {
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
-template<const int dim, typename BdryV_Scalar, typename BdryP, typename BdryPhi, typename BdryMu>
-class Bdry_2f0s
-{
-public:
- BdryV_Scalar bdryV_Scalar;
- BdryP bdryP;
- BdryPhi bdryPhi;
- BdryMu bdryMu;
- typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
- BdryV bdryV;
- BdryNS bdryNS;
- BdryCH bdryCH;
-
- Bdry_2f0s() :
- bdryV_Scalar(),
- bdryP(),
- bdryPhi(),
- bdryMu(),
- bdryV(bdryV_Scalar),
- bdryNS(bdryV,bdryP),
- bdryCH(bdryPhi,bdryMu)
- {
- }
-};
-
-//_________________________________________________________________________________________________________________________
-
template<typename GV, typename RF, typename EST_LOP, typename CH_GFS>
class RefinementIndicator
{
@@ -328,8 +141,8 @@ public:
using Z0 = Dune::PDELab::Backend::Vector<P0GFS,RF>;
Z0 z0;
- //typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
- //Z0DGF z0dgf;
+ typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
+ Z0DGF z0dgf;
//Estimator Operator
EST_LOP estlop;
@@ -347,7 +160,7 @@ public:
p0fem(Dune::GeometryTypes::simplex(GV::dimension)),
p0Gfs(gv,p0fem),
z0(p0Gfs,0.0),
- //z0dgf(p0Gfs,z0),
+ z0dgf(p0Gfs,z0),
estlop(),
estgo(chGfs,p0Gfs,estlop,MBE(1))
{
@@ -358,9 +171,15 @@ public:
}
template<typename CH_V>
- const Z0& calculateIndicator(const CH_V& chV)
+ const Z0& calculateIndicator(const GV& gv, const CH_V& chV, bool writeOutput)
{
estgo.residual(chV,z0);
+ if(writeOutput)
+ {
+ auto vtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+ vtkwriter->addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Z0DGF> >(z0dgf,"indicator"));
+ vtkwriter->write("RefinementIndicator",Dune::VTK::ascii);
+ }
return z0;
}
@@ -368,6 +187,6 @@ public:
void markGrid(Grid &grid, int verbosity=0)
{
if(verbosity>0) std::cout<<"mark elements"<<std::endl;
- Dune::PDELab::mark_grid(grid,z0,eta_refine,eta_coarsen, Rminlevel,Rmaxlevel,verbosity);
+ Dune::PDELab::mark_grid(grid,z0,eta_refine,eta_coarsen, Rminlevel,Rmaxlevel,verbosity,-1);
}
};
diff --git a/dune/phasefield/debug/vector_debug.hh b/dune/phasefield/debug/vector_debug.hh
index 1a81695..3c2a097 100644
--- a/dune/phasefield/debug/vector_debug.hh
+++ b/dune/phasefield/debug/vector_debug.hh
@@ -1,6 +1,6 @@
#pragma once
-#include "../myincludes.hh"
+#include <iostream>
#include <string>
diff --git a/dune/phasefield/ff_chns.hh b/dune/phasefield/ff_chns.hh
new file mode 100644
index 0000000..a8414da
--- /dev/null
+++ b/dune/phasefield/ff_chns.hh
@@ -0,0 +1,181 @@
+#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 P_FEM, typename PHI_FEM>
+class GFS_ff_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 conp, conphi, conmu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
+ VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs;
+
+ 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, V_GFS, P_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_ff_chns(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
+ conp(),conphi(),conmu(),
+ vGfs(es,vFem),
+ pGfs(es,pFem, conp),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs),
+ nsCC(),
+ chCC()
+ {
+ vGfs.name("velocity");
+ pGfs.name("pressure");
+ 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_ff_chns
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
+ PDGF pDgf;
+
+ 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_ff_chns(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
+ vSubGfs(gfs.ns),
+ vDgf(vSubGfs,nsV),
+ pSubGfs(gfs.ns),
+ pDgf(pSubGfs,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<VDGF> >(vDgf,"v"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu>
+class Ini_ff_chns
+{
+public:
+ IniV iniV;
+ IniP iniP;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ template<typename GV, typename Params>
+ Ini_ff_chns(const GV& gv, const Params& params) :
+ iniV(gv, params),
+ iniP(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniNS(iniV,iniP),
+ iniCH(iniPhi,iniMu)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<const int dim, typename BdryV_Scalar, typename BdryP, typename BdryPhi, typename BdryMu>
+class Bdry_ff_chns
+{
+public:
+ BdryV_Scalar bdryV_Scalar;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ Bdry_ff_chns() :
+ bdryV_Scalar(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryV(bdryV_Scalar),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu)
+ {
+ }
+};
diff --git a/dune/phasefield/ffs_chns_r.hh b/dune/phasefield/ffs_chns_r.hh
new file mode 100644
index 0000000..87ca6ea
--- /dev/null
+++ b/dune/phasefield/ffs_chns_r.hh
@@ -0,0 +1,232 @@
+#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>
+#include<dune/phasefield/utility/pf_3p_color_dgf.hh>
+
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM, typename U_FEM>
+class GFS_ffs_chns_r
+{
+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 conp, conphi, conmu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
+ VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs;
+
+ 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::GridFunctionSpace<ES, U_FEM, ConstraintsAssembler, VectorBackend> U_GFS;
+ U_GFS uGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS, U_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_ffs_chns_r(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem, const U_FEM& uFem) :
+ conp(),conphi(),conmu(),
+ vGfs(es,vFem),
+ pGfs(es,pFem, conp),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ uGfs(es,uFem, conmu),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs, phiGfs, muGfs, uGfs),
+ nsCC(),
+ chCC()
+ {
+ vGfs.name("velocity");
+ pGfs.name("pressure");
+ phiGfs.name("phasefield");
+ muGfs.name("potential");
+ uGfs.name("concentration");
+ }
+
+ 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_ffs_chns_r
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
+ PDGF pDgf;
+
+ 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;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<2> > Phi2SubGFS;
+ Phi2SubGFS phi2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Phi2SubGFS,CH_V> Phi2DGF;
+ Phi2DGF phi2Dgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<3> > Mu2SubGFS;
+ Mu2SubGFS mu2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Mu2SubGFS,CH_V> Mu2DGF;
+ Mu2DGF mu2Dgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<4> > USubGFS;
+ USubGFS uSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<USubGFS,CH_V> UDGF;
+ UDGF uDgf;
+
+ typedef PF_3P_Color_DGF<typename GFS::CH,CH_V,3> ColorDGF;
+ ColorDGF colorDgf;
+
+
+ DGF_ffs_chns_r(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
+ vSubGfs(gfs.ns),
+ vDgf(vSubGfs,nsV),
+ pSubGfs(gfs.ns),
+ pDgf(pSubGfs,nsV),
+
+ phiSubGfs(gfs.ch),
+ phiDgf(phiSubGfs,chV),
+ muSubGfs(gfs.ch),
+ muDgf(muSubGfs,chV),
+ phi2SubGfs(gfs.ch),
+ phi2Dgf(phi2SubGfs,chV),
+ mu2SubGfs(gfs.ch),
+ mu2Dgf(mu2SubGfs,chV),
+
+ uSubGfs(gfs.ch),
+ uDgf(uSubGfs,chV),
+
+ colorDgf(gfs.ch,chV)
+ {
+ }
+
+ template<typename VTKWRITER>
+ void addToVTKWriter(VTKWRITER& vtkwriter)
+ {
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Phi2DGF> >(phi2Dgf,"phi2"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Mu2DGF> >(mu2Dgf,"mu2"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<UDGF> >(uDgf,"u"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ColorDGF> >(colorDgf,"_coloring_"));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu, typename IniPhi2, typename IniMu2, typename IniU>
+class Ini_ffs_chns_r
+{
+public:
+ IniV iniV;
+ IniP iniP;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ IniPhi2 iniPhi2;
+ IniMu2 iniMu2;
+ IniU iniU;
+ typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu,IniPhi2,IniMu2,IniU> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ template<typename GV, typename Params>
+ Ini_ffs_chns_r(const GV& gv, const Params& params) :
+ iniV(gv, params),
+ iniP(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniPhi2(gv, params),
+ iniMu2(gv, params),
+ iniU(gv, params),
+ iniNS(iniV,iniP),
+ iniCH(iniPhi,iniMu,iniPhi2,iniMu2,iniU)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<const int dim, typename BdryV_Scalar, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2, typename BdryU>
+class Bdry_ffs_chns_r
+{
+public:
+ BdryV_Scalar bdryV_Scalar;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ BdryPhi2 bdryPhi2;
+ BdryMu2 bdryMu2;
+ BdryU bdryU;
+ typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu, BdryPhi2, BdryMu2, BdryU> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ Bdry_ffs_chns_r() :
+ bdryV_Scalar(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryPhi2(),
+ bdryMu2(),
+ bdryU(),
+ bdryV(bdryV_Scalar),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryPhi2,bdryMu2,bdryU)
+ {
+ }
+};
diff --git a/dune/phasefield/fs_chns_r.hh b/dune/phasefield/fs_chns_r.hh
new file mode 100644
index 0000000..f88ad39
--- /dev/null
+++ b/dune/phasefield/fs_chns_r.hh
@@ -0,0 +1,198 @@
+#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 P_FEM, typename PHI_FEM, typename U_FEM>
+class GFS_fs_chns_r
+{
+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 conp, conphi, conmu, conu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
+ VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs;
+
+ 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::GridFunctionSpace<ES, U_FEM, ConstraintsAssembler, VectorBackend> U_GFS;
+ U_GFS uGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS, U_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_fs_chns_r(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem, const U_FEM& uFem) :
+ conp(),conphi(),conmu(),
+ vGfs(es,vFem),
+ pGfs(es,pFem, conp),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ uGfs(es,uFem,conu),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs, uGfs),
+ nsCC(),
+ chCC()
+ {
+ vGfs.name("velocity");
+ pGfs.name("pressure");
+ phiGfs.name("phasefield");
+ muGfs.name("potential");
+ uGfs.name("concentration");
+ }
+
+ 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_fs_chns_r
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
+ PDGF pDgf;
+
+ 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;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<2> > USubGFS;
+ USubGFS uSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<USubGFS,CH_V> UDGF;
+ UDGF uDgf;
+
+
+ DGF_fs_chns_r(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
+ vSubGfs(gfs.ns),
+ vDgf(vSubGfs,nsV),
+ pSubGfs(gfs.ns),
+ pDgf(pSubGfs,nsV),
+ phiSubGfs(gfs.ch),
+ phiDgf(phiSubGfs,chV),
+ muSubGfs(gfs.ch),
+ muDgf(muSubGfs,chV),
+ uSubGfs(gfs.ch),
+ uDgf(uSubGfs,chV)
+ {
+ }
+
+ template<typename VTKWRITER>
+ void addToVTKWriter(VTKWRITER& vtkwriter)
+ {
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<UDGF> >(uDgf,"u"));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu, typename IniU>
+class Ini_fs_chns_r
+{
+public:
+ IniV iniV;
+ IniP iniP;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ IniU iniU;
+ typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu,IniU> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ template<typename GV, typename Params>
+ Ini_fs_chns_r(const GV& gv, const Params& params) :
+ iniV(gv, params),
+ iniP(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniU(gv, params),
+ iniNS(iniV,iniP),
+ iniCH(iniPhi,iniMu,iniU)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<const int dim, typename BdryV_Scalar, typename BdryP, typename BdryPhi, typename BdryMu, typename BdryU>
+class Bdry_fs_chns_r
+{
+public:
+ BdryV_Scalar bdryV_Scalar;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ BdryU bdryU;
+ typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> 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() :
+ bdryV_Scalar(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryU(),
+ bdryV(bdryV_Scalar),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryU)
+ {
+ }
+};
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index 52fb625..265d77f 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -1,7 +1,8 @@
#pragma once
-#include "myincludes.hh"
#include <cmath>
+#include <dune/pdelab/common/function.hh>
+#include <dune/phasefield/base_parameters.hh>
template<typename GV, typename RF, std::size_t dim_range, typename Params>
@@ -31,7 +32,35 @@ public:
ZeroInitial(const GV & gv, const Params & params) : ZeroInitialVector<GV,RF,1,Params>(gv,params) {}
};
-//_________________________________________________________________________________________________________________________
+//______________________________________________________________________________________________________________________
+
+#define CONST_INITIAL_U 0
+
+template<typename GV, typename RF, typename Params, const int type>
+class ConstInitial :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, ConstInitial<GV,RF,Params, type> >
+{
+public:
+ RF myval;
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ ConstInitial(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,ConstInitial<GV,RF,Params,type> > (gv)
+ {
+ if(type == CONST_INITIAL_U)
+ myval = params.initial.ptree.get("uConst",(RF)1);
+ else
+ myval = 0;
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ y = myval;
+ }
+};
+
+//______________________________________________________________________________________________________________________
template<typename GV, typename RF, typename Params>
class PhiInitial_RayleighTaylor :
@@ -55,4 +84,131 @@ public:
}
};
-//_________________________________________________________________________________________________________________________
+//______________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename Params, const int Phase>
+class PhiInitial_3pContactPoint :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_3pContactPoint<GV,RF,Params,Phase> >
+{
+ RF eps;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_3pContactPoint(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_3pContactPoint<GV,RF,Params,Phase> > (gv)
+ {
+ eps = params.eps;
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ y = 0.5*(1-tanh(1/eps/sqrt(2)*(x[1]-0.75-0.01*cos(x[0]*2*3.1415))));
+ y = 0.0;
+
+ RF d1 = x[0]-0.5;
+ RF yx = 0.5*(1+tanh(1/eps/sqrt(2)*d1));
+
+ RF d2 = x[1]-0.75;
+ if(Phase == 0)
+ d2 = -d2;
+ y[0]=0.5*(1+tanh(1/eps/sqrt(2)*d2))*yx;
+ }
+};
+
+//______________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename Params>
+class PhiInitial_Sphere :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Sphere<GV,RF,Params> >
+{
+ RF eps;
+ RF radius;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_Sphere(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Sphere<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 dist = sqrt( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5) );
+ y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-radius)));
+ }
+};
+
+//______________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename Params, int dim>
+class VelocityInitial_Inflow :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
+ VelocityInitial_Inflow<GV,RF,Params,dim> >
+{
+ const Params_Time& time;
+ RF meanflow;
+
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
+
+ typedef typename Traits::DomainType DomainType;
+ typedef typename Traits::RangeType RangeType;
+
+ VelocityInitial_Inflow(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,VelocityInitial_Inflow<GV,RF,Params,dim> > (gv)
+ , time(params.time)
+ {
+ meanflow = params.initial.ptree.get("MeanInflow",(RF)0.1);
+ }
+
+ template <typename T>
+ void setTime(T t){
+ time = t;
+ }
+
+ inline void evaluateGlobal(const DomainType & x, RangeType & y) const
+ {
+ y = 0.0;
+ if(x[0]<1e-7) //Inflow boundary
+ {
+ y[0] = std::max((RF)0,3/2 * meanflow * (x[1]+1)*(1-x[1]));
+ }
+ }
+};
+
+
+template<typename GV, typename RF, typename Params>
+class UInitial_Inflow :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, UInitial_Inflow<GV,RF,Params> >
+{
+ RF uConst;
+ RF uInflow;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ UInitial_Inflow(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,UInitial_Inflow<GV,RF,Params> > (gv)
+ {
+ uConst = params.initial.ptree.get("uConst",(RF)1.0);
+ uInflow = params.initial.ptree.get("uInflow",(RF)1.5);
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ y = uConst;
+ if(x[0]<1e-7)
+ {
+ y=uInflow;
+ }
+ }
+};
diff --git a/dune/phasefield/localoperator/fem_base_chns.hh b/dune/phasefield/localoperator/fem_base_chns.hh
new file mode 100644
index 0000000..a6ef66e
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_base_chns.hh
@@ -0,0 +1,298 @@
+#pragma once
+
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+
+template<typename GFS, typename Container, typename Entity, typename RangeFieldType>
+class GFS_LocalViewer
+{
+ typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
+ mutable LFS lfs;
+ typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
+ mutable LFSCache lfsCache;
+ typedef typename Container::template ConstLocalView<LFSCache> LView;
+ mutable LView lView;
+
+ mutable std::vector<RangeFieldType> lv;
+
+ 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)
+ {
+ if(entity == newEntity)
+ {
+ return false;
+ }
+ entity = newEntity;
+
+ lfs.bind(entity);
+ lfsCache.update();
+
+ lView.bind(lfsCache);
+ lView.read_sub_container(femspace,lv);
+ lView.unbind();
+
+ return true;
+ }
+
+ bool DoEvaluation(const Entity& newEntity)
+ {
+ if(entity == newEntity)
+ {
+ return false;
+ }
+ entity = newEntity;
+
+ lfs.bind(entity);
+ lfsCache.update();
+
+ lView.bind(lfsCache);
+ lView.read(lv);
+ lView.unbind();
+
+ return true;
+ }
+
+};
+
+
+// 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/chns_fem_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
similarity index 85%
rename from dune/phasefield/localoperator/chns_fem_cahnhilliard.hh
rename to dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
index 8318d15..ab36cc7 100644
--- a/dune/phasefield/localoperator/chns_fem_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh
@@ -1,18 +1,20 @@
#pragma once
-//#include "../myincludes.hh"
-
#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>
//Parameter class (for functions), FEM space and RangeField
template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
-class CHNS_FEM_CahnHilliard :
+class FEM_FF_CHNS_CahnHilliard :
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
- public Dune::PDELab::NumericalJacobianVolume<CHNS_FEM_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<CHNS_FEM_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FF_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FF_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
{
private:
typedef typename CHContainer::field_type RF;
@@ -55,7 +57,7 @@ public:
enum {doAlphaVolume=true};
//constructor
- CHNS_FEM_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsVOld_) :
+ FEM_FF_CHNS_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())
@@ -147,14 +149,14 @@ public:
{
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
r.accumulate(pfspace,i,factor*(
- (pf-pfOld)*phi[i]/param.dt
- +param.Mob(pfOld)*param.eps* (gradxi*gradphi[i][0])
+ (pf-pfOld)*phi[i]/param.time.dt
+ +param.mobility.eval(pfOld)*param.eps* (gradxi*gradphi[i][0])
-pf*(v*gradphi[i][0])
));
r.accumulate(xispace,i,factor*(
xi*phi[i]
- - 1/param.eps* (param.dWConvex(pf)+param.dWConcarve(pfOld))*phi[i]
+ - 1/param.eps* (param.dw.eval_dConvex(pf)+param.dw.eval_dConcarve(pfOld))*phi[i]
- param.eps/2 * (gradpf*gradphi[i][0])
));
}
diff --git a/dune/phasefield/localoperator/chns_fem_navierstokes.hh b/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
similarity index 87%
rename from dune/phasefield/localoperator/chns_fem_navierstokes.hh
rename to dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
index f8c8997..9461fda 100644
--- a/dune/phasefield/localoperator/chns_fem_navierstokes.hh
+++ b/dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh
@@ -1,18 +1,22 @@
#pragma once
-//#include "../myincludes.hh"
-
#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>
template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
- class CHNS_FEM_NavierStokes :
+ class FEM_FF_CHNS_NavierStokes :
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
- public Dune::PDELab::NumericalJacobianVolume<CHNS_FEM_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
- public Dune::PDELab::NumericalJacobianApplyVolume<CHNS_FEM_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FF_CHNS_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FF_CHNS_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
{
private:
typedef typename CHContainer::field_type RF;
@@ -69,7 +73,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
//TODO enum { doLambdaBoundary = true };
- CHNS_FEM_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
+ FEM_FF_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
NSGFS& nsGfs_, NSContainer& nsVOld_) :
param(param_),
chLfs(chGfs_), chLfsCache(chLfs), chView(chV_), chOldView(chVOld_),
@@ -179,10 +183,10 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
//Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
//g[1]= 20*xg[0]*(1-xg[0]);
g[1] = -5;
- RF rho=param.rho(pf);
- RF rhoOld=param.rho(pfOld);
+ RF rho=param.phys.rho_f(pf,1-pf);
+ RF rhoOld=param.phys.rho_f(pfOld,1-pfOld);
RF realPf = pf;
- pf=0.01+0.99*pf; //!!!!
+ //pf=0.01+0.99*pf; //!!!!
auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
for (size_t i=0; i<femspace.size(); i++)
@@ -206,10 +210,10 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
const auto v_nabla_v = v * jacv[k];
r.accumulate(child(vspace,k),i, factor*(
- ((rho+rhoOld)/2*v[k]-rhoOld*vOld[k])*vphi[i]/param.dt //TimeEvolution
- + ( rho* v_nabla_v + param.Mob(realPf)*gradxi[k]*(param.rho_1-param.rho_2) )* vphi[i] // Momentum transport
+ ((rho+rhoOld)/2*v[k]-rhoOld*vOld[k])*vphi[i]/param.time.dt //TimeEvolution
+ + ( rho* v_nabla_v + param.mobility.eval(realPf)*gradxi[k]*(param.phys.rho_1-param.phys.rho_2) )* vphi[i] // Momentum transport
- p * gradvphi[i][0][k] //Pressure
- + param.mu * (jacv[k]*gradvphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.phys.mu * (jacv[k]*gradvphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
- rho*g[k]*vphi[i] //Gravitation
+ param.SurfaceTensionScaling*pfOld*gradxi[k]*vphi[i] //Surface Tension
));
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh b/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
new file mode 100644
index 0000000..f8a5ecd
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
@@ -0,0 +1,239 @@
+#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>
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+ class FEM_FFS_CHNS_NavierStokes :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FFS_CHNS_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FFS_CHNS_NavierStokes<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;
+
+ 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:
+
+
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ FEM_FFS_CHNS_NavierStokes (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 = 2; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+ //const auto& myordering = lfsv.gridFunctionSpace().ordering();
+
+ const auto& pf1space = child(chLfs,_0);
+ const auto& xi1space = child(chLfs,_1);
+ const auto& pf2space = child(chLfs,_2);
+ const auto& xi2space = child(chLfs,_3);
+ 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();
+
+ chOldView.bind(chLfsCache);
+ chOldView.read(chOldLocal);
+ chOldView.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());
+
+ // 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::FieldVector<RF,dim> vOld(0.0);
+ Dune::FieldMatrix<RF,dim,dim> jacv(0.0);
+ RF p=0.0;
+ RF pf1=0.0;
+ RF pf2=0.0;
+ RF xi1=0.0;
+ RF xi2=0.0;
+ Dune::FieldVector<RF,dim> gradpf1(0.0);
+ Dune::FieldVector<RF,dim> gradpf2(0.0);
+ Dune::FieldVector<RF,dim> gradxi1(0.0);
+ Dune::FieldVector<RF,dim> gradxi2(0.0);
+
+ for (size_t i=0; i<femspace.size(); i++)
+ {
+
+ p += x(pspace,i)*phi[i];
+
+ pf1 += chLocal[pf1space.localIndex(i)]*phi[i];
+ gradpf1.axpy(chLocal[pf1space.localIndex(i)],gradphi[i][0]);
+ pf2 += chLocal[pf2space.localIndex(i)]*phi[i];
+ gradpf2.axpy(chLocal[pf2space.localIndex(i)],gradphi[i][0]);
+
+ xi1 += chLocal[xi1space.localIndex(i)]*phi[i];
+ gradxi1.axpy(chLocal[xi1space.localIndex(i)],gradphi[i][0]);
+ xi2 += chLocal[xi2space.localIndex(i)]*phi[i];
+ gradxi2.axpy(chLocal[xi2space.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]);
+ }
+ }
+ //std::cout << "v:" << v << std::endl;
+ //std::cout << "vOld:" << vOld << std::endl;
+ //std::cout << "jacv:" << jacv << std::endl;
+ //std::cout << "p:" << p << std::endl;
+
+ RF pf3=1-pf1-pf2;
+ RF pf_f=pf1+pf2+2*param.delta*pf3;
+ Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ gradpf_f.axpy((1-2*param.delta),gradpf1+gradpf2);
+
+ RF rho_f=param.phys.rho_f(pf1,pf2);
+ RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
+
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+ for (size_t i=0; i<femspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: (div v)*phi
+ //std::cout << -1.0* jacv[k][k]*phi[i]*factor << std::endl;
+ r.accumulate(pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);
+ }
+ }
+ for (size_t i=0; i<vfemspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // 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*(gradxi1[k]/param.tw.Sigma1 + gradxi2[k]/param.tw.Sigma2);
+ RF S = 0;
+ //S += xi1*(gradpf1[k]*pf_f - pf1 * gradpf_f[k])/pf_f;
+ //S += xi1*(gradpf2[k]*pf_f - pf2 * gradpf_f[k])/pf_f;
+ //S += 2*param.delta*(gradxi3k-gradxi2[k]-gradxi1[k]);
+ S += pf1*gradxi1[k] + pf2*gradxi2[k];
+ S = S*param.SurfaceTensionScaling*vphi[i];
+
+ 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
+ + param.phys.mu * (jacv[k]*gradvphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.vDis.eval(pf_f)*v[k]*vphi[i] //Velocity Dissipation in Solid
+ //- rho*g[k]*vphi[i] //Gravitation
+ + S //Surface Tension
+ )*factor);
+
+ }
+ }
+ }
+ }
+
+ };
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
new file mode 100644
index 0000000..6b4a53f
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
@@ -0,0 +1,226 @@
+#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>
+
+//Parameter class (for functions), FEM space and RangeField
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+class FEM_FFS_CHNS_R_CahnHilliard :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FFS_CHNS_R_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FFS_CHNS_R_CahnHilliard<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;
+
+ Param& param; // parameter functions
+
+ //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;
+
+public:
+ enum {doPatternVolume=true};
+ //enum {doLambdaVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_FFS_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())
+ {
+ }
+
+
+ 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& pf2space = child(lfsv,_2);
+ const auto& xi2space = child(lfsv,_3);
+ const auto& uspace = child(lfsv,_4);
+
+ 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();
+ }
+
+ 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());
+
+ // transform gradients of shape functions to real element
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto gradphi = makeJacobianContainer(femspace);
+ for (size_t i=0; i<femspace.size(); i++)
+ S.mv(gradphihat[i][0],gradphi[i][0]);
+
+ //compute current values of the unknowns
+ Dune::FieldVector<RF,dim> v(0.0);
+ RF pf1=0.0;
+ RF pf1Old=0.0;
+ RF xi1=0.0;
+
+ RF pf2=0.0;
+ RF pf2Old=0.0;
+ RF xi2=0.0;
+
+ RF u=0.0;
+ RF uOld=0.0;
+
+ Dune::FieldVector<RF,dim> gradpf1(0.0);
+ Dune::FieldVector<RF,dim> gradxi1(0.0);
+ Dune::FieldVector<RF,dim> gradpf2(0.0);
+ Dune::FieldVector<RF,dim> gradxi2(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];
+
+ gradpf2.axpy(x(pf2space,i),gradphi[i][0]);
+ pf2 += x(pf2space,i)*phi[i];
+ gradxi2.axpy(x(xi2space,i),gradphi[i][0]);
+ xi2 += x(xi2space,i)*phi[i];
+
+ gradu.axpy(x(uspace,i),gradphi[i][0]);
+ u += x(uspace,i)*phi[i];
+
+ pf1Old += chOldLocal[pf1space.localIndex(i)]*phi[i];
+ pf2Old += chOldLocal[pf2space.localIndex(i)]*phi[i];
+ 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];
+ }
+ }
+ RF pf3 = 1-pf1-pf2;
+ RF xi3 = -param.tw.Sigma3*(xi1/param.tw.Sigma1 + xi2/param.tw.Sigma2);
+
+ RF pf_c = pf1;
+ RF pf_c_Old = pf1Old;
+ 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++)
+ {
+ RF R1 = -param.tw.eval_q(pf1,pf2,pf3)/param.eps*
+ (param.reaction.eval(u)+param.phys.curvatureAlpha*(xi1-xi3))*phi[i];
+ RF J1 = param.mobility.eval(pf1,pf2,pf3)*param.eps/param.tw.Sigma1 * (gradxi1*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
+ +J1
+ -pf1*(v*gradphi[i][0])
+ -R1
+ ));
+
+ r.accumulate(pf2space,i,factor*(
+ (pf2-pf2Old)*phi[i]/param.time.dt
+ +param.mobility.eval(pf1,pf2,pf3)*param.eps/param.tw.Sigma2 * (gradxi2*gradphi[i][0])
+ -pf2*(v*gradphi[i][0])
+ ));
+
+ r.accumulate(xi1space,i,factor*(
+ xi1*phi[i]
+ - 1/param.eps* param.tw.eval_dpf1(pf1,pf2,pf3) *phi[i]
+ - param.eps * param.tw.Sigma1* (gradpf1*gradphi[i][0])
+ ));
+
+ r.accumulate(xi2space,i,factor*(
+ xi2*phi[i]
+ - 1/param.eps* param.tw.eval_dpf2(pf1,pf2,pf3) *phi[i]
+ - param.eps * param.tw.Sigma2* (gradpf2*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*(
+
+ //(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)*(v*gradphi[i][0])
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
new file mode 100644
index 0000000..845a371
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
@@ -0,0 +1,217 @@
+#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>
+
+//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> >
+{
+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;
+
+ Param& param; // parameter functions
+
+ //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;
+
+public:
+ enum {doPatternVolume=true};
+ //enum {doLambdaVolume=true};
+ enum {doAlphaVolume=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())
+ {
+ }
+
+
+ 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();
+ }
+
+ 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];
+
+ gradu.axpy(x(uspace,i),gradphi[i][0]);
+ u += x(uspace,i)*phi[i];
+
+ 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_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++)
+ {
+ 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]);
+
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(pf1space,i,factor*(
+ (pf1-pf1Old)*phi[i]/param.time.dt
+ +J1
+ //+(pf1*div_v + (gradpf1*v))*phi[i]
+ -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])
+ ));
+
+
+ // 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*(
+
+ //(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]
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
new file mode 100644
index 0000000..fe559da
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
@@ -0,0 +1,246 @@
+#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_R_InstatStokes :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_R_InstatStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_R_InstatStokes<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_FS_CHNS_R_InstatStokes (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/pf_3p_diff_estimator.hh b/dune/phasefield/localoperator/pf_3p_diff_estimator.hh
new file mode 100644
index 0000000..a2db287
--- /dev/null
+++ b/dune/phasefield/localoperator/pf_3p_diff_estimator.hh
@@ -0,0 +1,84 @@
+#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>
+
+template<typename CHGFS>
+class PF_3P_DIFF_ESTIMATOR :
+ public Dune::PDELab::LocalOperatorDefaultFlags
+{
+private:
+ typedef typename CHGFS::template Child<0>::Type FEM;
+ typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
+ Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+ enum {dim=LocalBasis::Traits::dimDomain};
+
+public:
+ enum {doPatternVolume=false};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ PF_3P_DIFF_ESTIMATOR()
+ {
+ }
+
+
+ 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
+ {
+ using namespace Dune::Indices;
+ const auto& pf1space = child(lfsu,_0);
+ const auto& pf2space = child(lfsu,_2);
+ typedef decltype(makeZeroBasisFieldValue(pf1space)) RF;
+
+ auto geo = eg.geometry();
+ const auto& rElement = Dune::ReferenceElements<double, dim>::general(geo.type());
+ RF u1max=0.0;
+ RF u1min=0.0;
+ RF u2max=0.0;
+ RF u2min=0.0;
+ RF u3max=0.0;
+ RF u3min=0.0;
+ // Dune::FieldVector<RF, dim> xg = eg.geometry().global( rElement.position(0, dim) );
+ for (int i=0; i<geo.corners(); i++)
+ {
+ RF u1=0.0;
+ RF u2=0.0;
+ auto& phihat = cache.evaluateFunction(rElement.position(i, dim),pf1space.finiteElement().localBasis());
+ for (size_t j=0; j<pf1space.size(); j++)
+ {
+ u1 += x(pf1space,j)*phihat[j];
+ u2 += x(pf2space,j)*phihat[j];
+ }
+ RF u3=1-u1-u2;
+ if(i==0)
+ {
+ u1max = u1;
+ u1min = u1;
+ u2max = u2;
+ u2min = u2;
+ u3max = u3;
+ u3min = u3;
+ }
+ else
+ {
+ if(u1>u1max) u1max = u1;
+ if(u1<u1min) u1min = u1;
+ if(u2>u2max) u2max = u2;
+ if(u2<u2min) u2min = u2;
+ if(u3>u3max) u3max = u3;
+ if(u3<u3min) u3min = u3;
+ }
+ //xg = eg.geometry().global( rElement.position(i, dim) );
+ }
+ r.accumulate(lfsv,0,0.5*((u1max-u1min)+(u2max-u2min)+(u3max-u3min)));//(xg[0]>0.25&&xg[0]<0.75)? (umax-umin):0 );
+ }
+
+};
diff --git a/src/test_2f0s/acfs_fem_ac_estimator.hh b/dune/phasefield/localoperator/pf_diff_estimator.hh
similarity index 82%
rename from src/test_2f0s/acfs_fem_ac_estimator.hh
rename to dune/phasefield/localoperator/pf_diff_estimator.hh
index af958ce..b2afcea 100644
--- a/src/test_2f0s/acfs_fem_ac_estimator.hh
+++ b/dune/phasefield/localoperator/pf_diff_estimator.hh
@@ -1,13 +1,17 @@
#pragma once
-#include "myincludes.hh"
#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>
template<typename CHGFS>
-class ACFS_FEM_AC_Estimator :
+class PF_DIFF_ESTIMATOR :
public Dune::PDELab::LocalOperatorDefaultFlags
{
private:
@@ -18,11 +22,10 @@ private:
public:
enum {doPatternVolume=false};
- //enum {doLambdaVolume=true};
enum {doAlphaVolume=true};
//constructor
- ACFS_FEM_AC_Estimator()
+ PF_DIFF_ESTIMATOR()
{
}
diff --git a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
index bd1ac1f..4b79d18 100644
--- a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
+++ b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
@@ -1,17 +1,7 @@
#pragma once
-#include "../myincludes.hh"
-#include <dune/istl/io.hh>
-#include <dune/istl/operators.hh>
-#include <dune/istl/owneroverlapcopy.hh>
-#include <dune/istl/paamg/amg.hh>
-#include <dune/istl/paamg/pinfo.hh>
-#include <dune/istl/preconditioners.hh>
-#include <dune/istl/scalarproducts.hh>
-#include <dune/istl/solvercategory.hh>
-#include <dune/istl/solvers.hh>
-#include <dune/istl/superlu.hh>
+#include<dune/pdelab/backend/istl.hh>
template<class GO,
template<class,class,class,int> class Preconditioner,
diff --git a/dune/phasefield/utility/pf_3p_color_dgf.hh b/dune/phasefield/utility/pf_3p_color_dgf.hh
new file mode 100644
index 0000000..f171697
--- /dev/null
+++ b/dune/phasefield/utility/pf_3p_color_dgf.hh
@@ -0,0 +1,112 @@
+#pragma once
+
+#include <cmath>
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+
+
+template<typename T, typename X, std::size_t dimR>
+class PF_3P_Color_DGF
+ : public Dune::PDELab::GridFunctionBase<
+ Dune::PDELab::GridFunctionTraits<
+ typename T::Traits::GridViewType,
+ typename T::template Child<0>::Type::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeFieldType,
+ dimR,
+ Dune::FieldVector<
+ typename T::template Child<0>::Type::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeFieldType,
+ dimR
+ >
+ >,
+ PF_3P_Color_DGF<T,X, dimR>
+ >
+{
+ typedef T GFS;
+
+ typedef Dune::PDELab::GridFunctionBase<
+ Dune::PDELab::GridFunctionTraits<
+ typename T::Traits::GridViewType,
+ typename T::template Child<0>::Type::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeFieldType,
+ dimR,
+ Dune::FieldVector<
+ typename T::template Child<0>::Type::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeFieldType,
+ dimR
+ >
+ >,
+ PF_3P_Color_DGF<T,X,dimR>
+ > BaseT;
+
+public:
+ typedef typename BaseT::Traits Traits;
+ typedef typename T::template Child<0>::Type ChildType;
+ typedef typename ChildType::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeFieldType RF;
+ typedef typename ChildType::Traits::FiniteElementType
+ ::Traits::LocalBasisType::Traits::RangeType RT;
+
+
+ PF_3P_Color_DGF(const GFS& gfs, const X& x_)
+ : pgfs(stackobject_to_shared_ptr(gfs))
+ , lfs(gfs)
+ , lfs_cache(lfs)
+ , x_view(x_)
+ , xl(gfs.maxLocalSize())
+ , yb(gfs.maxLocalSize())
+ { }
+
+ inline void evaluate (const typename Traits::ElementType& e,
+ const typename Traits::DomainType& x,
+ typename Traits::RangeType& y) const
+ {
+ lfs.bind(e);
+ lfs_cache.update();
+ x_view.bind(lfs_cache);
+ x_view.read(xl);
+ x_view.unbind();
+
+ using namespace Dune::Indices;
+ const auto& pf1space = child(lfs,_0);
+ const auto& pf2space = child(lfs,_2);
+
+ pf1space.finiteElement().localBasis().
+ evaluateFunction(x,yb);
+ RF p1 = 0;
+ RF p2 = 0;
+
+ for (unsigned int i=0; i<yb.size(); i++)
+ {
+ p1 += xl[pf1space.localIndex(i)]*yb[i];
+ p2 += xl[pf2space.localIndex(i)]*yb[i];
+ }
+ RF p3 = 1-p1-p2;
+ p1=fmax(0,fmin(1,p1));
+ p2=fmax(0,fmin(1,p2));
+ p3=fmax(0,fmin(1,p3));
+ RF pmax=fmax(p1, fmax(p2,p3));
+ y[0]=p1/pmax;
+ y[1]=p2/pmax;
+ y[2]=p3/pmax;
+ }
+
+ inline const typename Traits::GridViewType& getGridView () const
+ {
+ return pgfs->gridView();
+ }
+
+private:
+ typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
+ typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
+ typedef typename X::template ConstLocalView<LFSCache> XView;
+ std::shared_ptr<GFS const> pgfs;
+ mutable LFS lfs;
+ mutable LFSCache lfs_cache;
+ mutable XView x_view;
+ mutable std::vector<RF> xl;
+ mutable std::vector<RT> yb;
+};
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index bf71666..fb18a38 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,5 @@
add_executable("dune-phasefield" dune-phasefield.cc)
target_link_dune_default_libraries("dune-phasefield")
add_subdirectory(test_2f0s)
+add_subdirectory(test_3pchns_r)
+add_subdirectory(test_2pchns_r)
diff --git a/src/test_2f0s/chns.ini b/src/test_2f0s/2f0s.ini
similarity index 59%
rename from src/test_2f0s/chns.ini
rename to src/test_2f0s/2f0s.ini
index 504f780..2763223 100644
--- a/src/test_2f0s/chns.ini
+++ b/src/test_2f0s/2f0s.ini
@@ -1,29 +1,26 @@
-[physics]
+[Physics]
uAst = 2
D = 0.01
-rho = 1 #only used in 2p
+rho1 = 1.0
+rho2 = 1.8
mu=0.005
-fscaling=3
SurfaceTensionScaling=0.1
-[3p]
-Sigma1=1
-Sigma2=1
-Sigma3=10
-rho1=1.8
-rho2=1
[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_53k_per.msh
-
#filename = grids/square_2k_del_per.msh
+
+[Phasefield]
+eps=0.03
+Mobility=0.03
+
[Time]
tMax = 30
dt = 0.015
@@ -33,20 +30,22 @@ dt = 0.015
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
-[problem]
-eps=0.03
-Mobility=0.03
-MobilityReduced=0.001
-URegularization=0.01
-delta=0.03
-
-[initial]
-u=1
-vmax=0
-vmaxtime=5
-
[Refinement]
etaRefine = 0.15
-etaCoarsen = 0.05
+etaCoarsen = 0.14
maxRefineLevel = 2
-minRefineLevel = 0
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/test_2f0s/CMakeLists.txt b/src/test_2f0s/CMakeLists.txt
index fd193ff..e0ac545 100644
--- a/src/test_2f0s/CMakeLists.txt
+++ b/src/test_2f0s/CMakeLists.txt
@@ -1,2 +1,2 @@
add_executable(chns_equation chns_equation.cc)
-dune_symlink_to_source_files(FILES grids chns.ini)
+dune_symlink_to_source_files(FILES grids 2f0s.ini)
diff --git a/src/test_2f0s/chns_driver.hh b/src/test_2f0s/chns_driver.hh
index e193058..9c93f45 100644
--- a/src/test_2f0s/chns_driver.hh
+++ b/src/test_2f0s/chns_driver.hh
@@ -1,33 +1,17 @@
#pragma once
-#include "myincludes.hh"
-//#include <dune/xt/grid/grids.hh>
-//#include <dune/xt/grid/view/periodic.hh>
-#include <dune/pdelab/localoperator/stokesparameter.hh>
-#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
-#include <dune/pdelab/localoperator/taylorhoodnavierstokes.hh>
-#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
-#include <dune/pdelab/adaptivity/adaptivity.hh>
-
-//#include "acfs_fem_allencahn.hh"
-//#include "acfs_fem_cahnhilliard.hh"
-//#include "acfs_fem_navierstokes.hh"
-#include "acfs_fem_ac_estimator.hh"
-//#include "threephase_fem_cahnhilliard.hh"
-//#include "threephase_fem_navierstokes.hh"
-//#include "threephase_dgf_coloring.hh"
-#include <dune/phasefield/localoperator/chns_fem_cahnhilliard.hh>
-#include <dune/phasefield/localoperator/chns_fem_navierstokes.hh>
-
-#include <dune/phasefield/acfs_problem.hh>
-#include <dune/phasefield/acfs_ibvp.hh>
+
+#include <dune/phasefield/localoperator/fem_ff_chns_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_ff_chns_navierstokes.hh>
+#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
+
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/ff_chns.hh>
#include <dune/phasefield/initial_fn.hh>
#include <dune/phasefield/debug/vector_debug.hh>
-//template <class A>
-//using MySolver = Dune::RestartedGMResSolver<A, A, A>;
+#include <dune/phasefield/base_parameters.hh>
//_________________________________________________________________________________________________________________________
@@ -39,27 +23,18 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//std::bitset<2> periodic_directions(std::string("01"));
//typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
GV gv = grid.leafGridView();
- // dimension and important types
const int dim = GV::dimension;
- typedef double RF; // type for computations
-
-
-
- RF initialvmax = ptree.get("initial.vmax",(RF)1);
- RF vmaxtime = ptree.get("initial.vmaxtime",(RF)1);
-
-
- TimestepManager timestepManager(ptree.sub("Time"),2);
-
- Problem<RF> problem(ptree.get("temporal.tau",(RF)0.1),ptree);
- problem.setTime(timestepManager.t);
+ typedef double RF;
+ typedef Params_ff<RF, DoubleWell_poly<RF>, Mobility_Const<RF>> Parameters;
+ Parameters param(ptree);
+ TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
//_________________________________________________________________________________________________________________________
// Make grid function spaces
- typedef GFS_2f0s<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+ typedef GFS_ff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
GFS gfs(es,vFem,pFem,phiFem);
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
@@ -68,7 +43,7 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
//Constraints
- typedef Bdry_2f0s<dim, Dune::PDELab::DirichletConstraintsParameters,
+ typedef Bdry_ff_chns<dim, Dune::PDELab::DirichletConstraintsParameters,
Dune::PDELab::NoDirichletConstraintsParameters,
Dune::PDELab::NoDirichletConstraintsParameters,
Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
@@ -88,22 +63,22 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
//Initial Values
- typedef Ini_2f0s< ZeroInitialVector<GV,RF,dim,Problem<RF>>,
- ZeroInitial<GV,RF,Problem<RF>>,
- PhiInitial_RayleighTaylor<GV,RF,Problem<RF>>,
- ZeroInitial<GV,RF,Problem<RF>> > Initial;
- Initial initial(gv,problem);
+ typedef Ini_ff_chns< ZeroInitialVector<GV,RF,dim,Parameters>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_RayleighTaylor<GV,RF,Parameters>,
+ ZeroInitial<GV,RF,Parameters> > Initial;
+ Initial initial(gv,param);
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
//_________________________________________________________________________________________________________________________
// LocalOperators
- typedef CHNS_FEM_NavierStokes<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
- NS_LOP nsLop(problem, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+ typedef FEM_FF_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+ NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
- typedef CHNS_FEM_CahnHilliard<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
- CH_LOP chLop(problem, gfs.ch, chVOld, gfs.ns, nsV);
+ typedef FEM_FF_CHNS_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+ CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
@@ -111,7 +86,7 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
// Construct discrete Grid Functions
- typedef DGF_2f0s<GFS, NS_V, CH_V> DGF;
+ typedef DGF_ff_chns<GFS, NS_V, CH_V> DGF;
DGF dgf(gfs, nsV, chV);
//_________________________________________________________________________________________________________________________
@@ -119,22 +94,21 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
std::string filename=ptree.get("output.filename","output");
struct stat st;
if( stat( filename.c_str(), &st ) != 0 )
- {
- int stat = 0;
- stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
- if( stat != 0 && stat != -1)
- std::cout << "Error: Cannot create directory "
- << filename << std::endl;
- }
+ {
+ int stat = 0;
+ stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << filename << std::endl;
+ }
auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
dgf.addToVTKWriter(vtkwriter);
- vtkwriter.write(timestepManager.t,Dune::VTK::ascii);
+ vtkwriter.write(param.time.t,Dune::VTK::ascii);
//_________________________________________________________________________________________________________________________
- printVector(chV,50,"chV");
- printVector(nsV,50,"nsV");
+ printVector(chV,10,"chV");
+ printVector(nsV,10,"nsV");
//____________________________________________________________________________________________________________________
@@ -143,8 +117,6 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
nsVOld = nsV;
while (!timestepManager.isFinished())
{
- problem.setTime(timestepManager.t);
-
// Grid Operators
typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
@@ -159,25 +131,14 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
-
+ // Nonlinear solver
typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
PDESolver1 ns_newton(nsGo,nsV,ls_ns);
- ns_newton.setReassembleThreshold(0.0);
- ns_newton.setVerbosityLevel(2);
- ns_newton.setMaxIterations(50);
- ns_newton.setLineSearchMaxIterations(30);
-
- ns_newton.setAbsoluteLimit(1e-9);
- ns_newton.setReduction(1e-7);
+ ns_newton.setParameters(ptree.sub("NS_Newton"));
typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
PDESolver2 ch_newton(chGo,chV,ls_ch);
- ch_newton.setReassembleThreshold(0.0);
- ch_newton.setVerbosityLevel(1);
- ch_newton.setMaxIterations(50);
- ch_newton.setLineSearchMaxIterations(50);
-
-
+ ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
//_________________________________________________________________________________________________________________________
// do time step
@@ -203,9 +164,9 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
// mark elements for refinement
- typedef RefinementIndicator<GV,RF,ACFS_FEM_AC_Estimator<CH_GFS>,CH_GFS> RIND;
+ typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
- rInd.calculateIndicator(chV);
+ rInd.calculateIndicator(gv,chV,false);
rInd.markGrid(grid,2);
printVector(chV,3,"chV");
@@ -213,19 +174,22 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
// do refinement
std::cout << "adapt grid and solution" << std::endl;
- auto chTransfer = transferSolutions(gfs.ch,2,chV,chVOld);
- auto nsTransfer = transferSolutions(gfs.ns,2,nsV,nsVOld);
+ auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+ auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
std::cout << "adapt grid and solution ended" << std::endl;
printVector(chV,3,"chV");
printVector(nsV,3,"nsV");
+ //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
// recompute constraints
gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
- vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+ //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
}
}
diff --git a/src/test_2f0s/chns_equation.cc b/src/test_2f0s/chns_equation.cc
index a058dbc..1a1b97d 100644
--- a/src/test_2f0s/chns_equation.cc
+++ b/src/test_2f0s/chns_equation.cc
@@ -37,12 +37,12 @@ int main(int argc, char** argv)
// open ini file
Dune::ParameterTree ptree;
Dune::ParameterTreeParser ptreeparser;
- ptreeparser.readINITree("chns.ini",ptree);
+ ptreeparser.readINITree("2f0s.ini",ptree);
ptreeparser.readOptions(argc,argv,ptree);
// read ini file
// const int dim = ptree.get<int>("grid.dim");
- const int refinement = ptree.get<int>("domain.level");
+// const int refinement = ptree.get<int>("domain.level");
if (true/*dim==2*/)
diff --git a/src/test_3pchns_r/3p_driver.hh b/src/test_3pchns_r/3p_driver.hh
new file mode 100644
index 0000000..97347a3
--- /dev/null
+++ b/src/test_3pchns_r/3p_driver.hh
@@ -0,0 +1,216 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include <dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh>
+#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/ffs_chns_r.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_3p (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+ typedef typename Grid::LeafGridView GV;
+ //std::bitset<2> periodic_directions(std::string("01"));
+ //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+ GV gv = grid.leafGridView();
+ const int dim = GV::dimension;
+ typedef double RF;
+
+
+ typedef Params_ffs_r< RF,
+ TripleWell<RF,DoubleWell_limited<RF>>,
+ Mobility_Const<RF>,
+ Reaction_Quadratic_Limited<RF>,
+ VelDissipation_Quadratic<RF> > Parameters;
+ Parameters param(ptree);
+ TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make grid function spaces
+
+
+ typedef GFS_ffs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+ GFS gfs(es,vFem,pFem,phiFem,phiFem);
+ typedef typename GFS::NS NS_GFS;
+ typedef typename GFS::CH CH_GFS;
+
+
+ //_________________________________________________________________________________________________________________________
+ //Constraints
+
+ typedef Bdry_ffs_chns_r<dim, Dune::PDELab::DirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
+ Bdry bdry;
+ gfs.updateConstraintsContainer(bdry);
+
+ //_________________________________________________________________________________________________________________________
+ //Coefficient vectors
+
+ using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
+ using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
+ CH_V chV(gfs.ch);
+ NS_V nsV(gfs.ns);
+ CH_V chVOld = chV;
+ NS_V nsVOld = nsV;
+
+ //_________________________________________________________________________________________________________________________
+ //Initial Values
+
+ typedef Ini_ffs_chns_r< ZeroInitialVector<GV,RF,dim,Parameters>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_3pContactPoint<GV,RF,Parameters,0>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_3pContactPoint<GV,RF,Parameters,1>,
+ ZeroInitial<GV,RF,Parameters>,
+ ZeroInitial<GV,RF,Parameters> > Initial;
+ Initial initial(gv,param);
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+ //_________________________________________________________________________________________________________________________
+ // LocalOperators
+
+ typedef FEM_FFS_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+ NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+
+ typedef FEM_FFS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+ CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+ MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+
+
+ //_________________________________________________________________________________________________________________________
+ // Construct discrete Grid Functions
+ typedef DGF_ffs_chns_r<GFS, NS_V, CH_V> DGF;
+ DGF dgf(gfs, nsV, chV);
+
+ //_________________________________________________________________________________________________________________________
+ // prepare VTK writer
+ std::string filename=ptree.get("output.filename","output");
+ struct stat st;
+ if( stat( filename.c_str(), &st ) != 0 )
+ {
+ int stat = 0;
+ stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << filename << std::endl;
+ }
+ auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+ Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+ dgf.addToVTKWriter(vtkwriter);
+
+ vtkwriter.write(param.time.t,Dune::VTK::ascii);
+
+//_________________________________________________________________________________________________________________________
+ printVector(chV,10,"chV");
+ printVector(nsV,10,"nsV");
+
+
+ //____________________________________________________________________________________________________________________
+ //int i_it=0;
+ chVOld = chV;
+ nsVOld = nsV;
+ while (!timestepManager.isFinished())
+ {
+
+ // Grid Operators
+ typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+ NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
+
+ typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+ CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
+
+ // Linear solver
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+ <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
+
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+ <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
+
+ // Nonlinear solver
+ typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
+ PDESolver1 ns_newton(nsGo,nsV,ls_ns);
+ ns_newton.setParameters(ptree.sub("NS_Newton"));
+
+ typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
+ PDESolver2 ch_newton(chGo,chV,ls_ch);
+ ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+ //_________________________________________________________________________________________________________________________
+ // do time step
+
+
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+ int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+ if(result <0)
+ {
+ return;
+ }
+ else if(result == 0)
+ {
+ nsV = nsVOld;
+ chV = chVOld;
+ continue;
+ }
+
+ chVOld = chV;
+ nsVOld = nsV;
+
+ //_________________________________________________________________________________________________________________________
+
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chV,false);
+ rInd.markGrid(grid,2);
+
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
+
+ // do refinement
+ std::cout << "adapt grid and solution" << std::endl;
+ auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+ auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
+ Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
+ std::cout << "adapt grid and solution ended" << std::endl;
+
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
+
+ //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
+ // recompute constraints
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
+
+ //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+ }
+}
diff --git a/src/test_3pchns_r/3p_equation.cc b/src/test_3pchns_r/3p_equation.cc
new file mode 100644
index 0000000..b17ed0b
--- /dev/null
+++ b/src/test_3pchns_r/3p_equation.cc
@@ -0,0 +1,116 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include<iostream>
+// dune includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+//#if HAVE_UG
+//#include<dune/grid/uggrid.hh>
+//#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/alugrid/dgf.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// pdelab includes
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+
+#include "3p_driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+
+
+
+int main(int argc, char** argv)
+{
+ #if !HAVE_SUPERLU
+ std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
+ exit(1);
+ #endif
+
+ #if HAVE_UG
+ try{
+ // Maybe initialize Mpi
+ Dune::MPIHelper&
+ helper = Dune::MPIHelper::instance(argc, argv);
+ if(Dune::MPIHelper::isFake)
+ std::cout<< "This is a sequential program." << std::endl;
+ else
+ std::cout << "Parallel code run on "
+ << helper.size() << " process(es)" << std::endl;
+
+ // open ini file
+ Dune::ParameterTree ptree;
+ Dune::ParameterTreeParser ptreeparser;
+ ptreeparser.readINITree("chns.ini",ptree);
+ ptreeparser.readOptions(argc,argv,ptree);
+
+ // read ini file
+// const int dim = ptree.get<int>("grid.dim");
+ const int refinement = ptree.get<int>("domain.level");
+
+
+ if (true/*dim==2*/)
+ {
+
+ typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+ //typedef Dune::UGGrid<2> Grid;
+ typedef Grid::ctype DF;
+
+ std::vector<int> boundary_index_map;
+ std::vector<int> element_index_map;
+
+ std::string filename = ptree.get("domain.filename",
+ "turbtube2d.msh");
+ Dune::GridFactory<Grid> factory;
+ if (helper.rank()==0)
+ Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
+ std::shared_ptr<Grid> gridp(factory.createGrid());
+ //gridp->globalRefine(refinement);
+ //gridp->loadBalance();
+ //std::cout << " after load balance /" << helper.rank() << "/ " << gridp->size(0) << std::endl;
+
+ //typedef Grid::LeafGridView LGV;
+
+ //std::bitset<2> periodic_directions(std::string("01"));
+ //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+ //GV gv(gridp->leafGridView(),periodic_directions);
+
+ typedef Grid::LeafGridView GV;
+ GV gv = gridp->leafGridView();
+
+ std::cout << "OK" << std::endl;
+
+ using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+ ES es(gv);
+
+ const int k = 2;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
+ const int l = 1;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+ V_FEM vFem(es);
+ P_FEM pFem(es);
+ PHI_FEM phiFem(es);
+
+ driver_3p(*gridp, es,vFem,pFem,phiFem, ptree, helper);
+ }
+
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 1;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ return 1;
+ }
+ #endif
+}
diff --git a/src/test_3pchns_r/CMakeLists.txt b/src/test_3pchns_r/CMakeLists.txt
new file mode 100644
index 0000000..7669afd
--- /dev/null
+++ b/src/test_3pchns_r/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(3p_equation 3p_equation.cc)
+dune_symlink_to_source_files(FILES grids chns.ini)
diff --git a/src/test_3pchns_r/chns.ini b/src/test_3pchns_r/chns.ini
new file mode 100644
index 0000000..998a2eb
--- /dev/null
+++ b/src/test_3pchns_r/chns.ini
@@ -0,0 +1,58 @@
+[physics]
+uAst = 2
+D = 0.01
+rho1 = 1.0
+rho2 = 1.0
+mu=0.005
+ReactionRate=3
+#ReactionLimiter = 1000
+SurfaceTensionScaling=0.1
+#VelDissipationRate = 10
+CurvatureAlpha=0.03
+
+[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_53k_per.msh
+#filename = grids/square_2k_del_per.msh
+
+[Phasefield]
+eps=0.03
+delta=0.03
+DoubleWellDelta=0.5
+Mobility=0.03
+Sigma1=1
+Sigma2=1
+Sigma3=10
+
+[Time]
+tMax = 30
+dt = 0.015
+#dtMax = 0.015
+dtMin = 0.0015
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Refinement]
+etaRefine = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 3
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/test_3pchns_r/myincludes.hh b/src/test_3pchns_r/myincludes.hh
new file mode 100644
index 0000000..86fe412
--- /dev/null
+++ b/src/test_3pchns_r/myincludes.hh
@@ -0,0 +1,59 @@
+#pragma once
+
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<math.h>
+#include<iostream>
+// dune-common includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+// dune-geometry includes
+#include<dune/geometry/referenceelements.hh>
+#include<dune/geometry/quadraturerules.hh>
+// dune-grid includes
+#include<dune/grid/onedgrid.hh>
+#include<dune/grid/yaspgrid.hh>
+#include<dune/grid/utility/structuredgridfactory.hh>
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/grid/io/file/gmshreader.hh>
+#if HAVE_UG
+#include<dune/grid/uggrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/alugrid/dgf.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// dune-istl included by pdelab
+// dune-pdelab includes
+#include<dune/pdelab/common/function.hh>
+#include<dune/pdelab/common/vtkexport.hh>
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/finiteelementmap/qkfem.hh>
+#include<dune/pdelab/finiteelementmap/p0fem.hh>
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/common/constraintsparameters.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+#include<dune/pdelab/function/callableadapter.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/subspace.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include<dune/pdelab/gridfunctionspace/interpolate.hh>
+#include<dune/pdelab/gridfunctionspace/vtk.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/gridoperator/onestep.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/backend/istl.hh>
+#include<dune/pdelab/stationary/linearproblem.hh>
+#include<dune/pdelab/instationary/onestep.hh>
+#include<dune/pdelab/newton/newton.hh>
+#include<dune/pdelab/adaptivity/adaptivity.hh>
+
+//#include <dune/xt/grid/grids.hh>
+//#include <dune/xt/grid/view/periodic.hh>
--
GitLab