diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index dc7244c6231485363e1a06369373fbb897fa9325..fbd474aaf9841439791be727412848e90d7195a4 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -253,8 +253,50 @@ void interpolateToBdry(const GFS& gfs, const NS_BDRY_FN& nsBdryFn, const CH_BDRY
Dune::PDELab::copy_constrained_dofs(gfs.chCC, chInterp, chV);
}
+template<typename GFS, typename CHNS_BDRY_FN, typename CHNS_V>
+void interpolateToBdry_monolithic(const GFS& gfs, const CHNS_BDRY_FN& chnsBdryFn, CHNS_V& chnsV)
+{
+ CHNS_V chnsInterp(gfs.chns);
+ Dune::PDELab::interpolate(chnsBdryFn, gfs.chns, chnsInterp);
+ Dune::PDELab::copy_constrained_dofs(gfs.chnsCC, chnsInterp, chnsV);
+}
+
//_________________________________________________________________________________________________________________________
+template<class Grid, class GV, class GFS1, class X1>
+void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, int int_order)
+{
+ typedef Dune::PDELab::L2Projection<GFS1,X1> Projection1;
+ Projection1 projection1(gfs1,int_order);
+
+ Dune::PDELab::GridAdaptor<Grid,GFS1,X1,Projection1> grid_adaptor1(gfs1);
+
+ // prepare the grid for refinement
+ grid.preAdapt();
+
+ // save solution
+ typename Dune::PDELab::GridAdaptor<Grid,GFS1,X1,Projection1>::MapType transferMap1;
+ grid_adaptor1.backupData(grid,gfs1,projection1,x1,transferMap1);
+
+ // adapt the grid
+
+ std::cout << "Grid vor adapt";
+ grid.adapt();
+ std::cout << "Grid nach adapt";
+ gv.update();
+ std::cout << "nach GV adapt";
+
+ // update the function spaces
+ gfs1.update(true);
+
+ // interpolate solution
+ x1 = X1(gfs1,0.0);
+ grid_adaptor1.replayData(grid,gfs1,projection1,x1,transferMap1);
+
+ // clean up
+ grid.postAdapt();
+}
+
template<class Grid, class GV, class GFS1, class GFS2, class X1, class X2>
void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, GFS2& gfs2, X2& x2, int int_order)
{
@@ -297,6 +339,42 @@ void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, GFS2& gfs2, X2& x2,
grid.postAdapt();
}
+template<class Grid, class GV, class GFS1, class X1>
+void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, X1& x1b, int int_order)
+{
+ typedef Dune::PDELab::L2Projection<GFS1,X1> Projection1;
+ Projection1 projection1(gfs1,int_order);
+
+ Dune::PDELab::GridAdaptor<Grid,GFS1,X1,Projection1> grid_adaptor1(gfs1);
+ Dune::PDELab::GridAdaptor<Grid,GFS1,X1,Projection1> grid_adaptor1b(gfs1);
+
+ // prepare the grid for refinement
+ grid.preAdapt();
+
+ // save solution
+ typename Dune::PDELab::GridAdaptor<Grid,GFS1,X1,Projection1>::MapType transferMap1;
+ grid_adaptor1.backupData(grid,gfs1,projection1,x1,transferMap1);
+ grid_adaptor1b.backupData(grid,gfs1,projection1,x1b,transferMap1);
+
+ // adapt the grid
+
+ std::cout << "Grid vor adapt";
+ grid.adapt();
+ std::cout << "Grid nach adapt";
+ gv.update();
+ std::cout << "nach GV adapt";
+
+ // update the function spaces
+ gfs1.update(true);
+
+ // interpolate solution
+ x1 = X1(gfs1,0.0);
+ grid_adaptor1.replayData(grid,gfs1,projection1,x1,transferMap1);
+ x1b = X1(gfs1,0.0);
+ grid_adaptor1b.replayData(grid,gfs1,projection1,x1b,transferMap1);
+ // clean up
+ grid.postAdapt();
+}
template<class Grid, class GV, class GFS1, class GFS2, class X1, class X2>
void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, X1& x1b, GFS2& gfs2, X2& x2, X2& x2b, int int_order)
@@ -491,6 +569,27 @@ public:
}
}
+ template<typename Grid, typename CHNS_V, typename CHNS_GFS>
+ void myRefineGrid(Grid &grid, GV& gv, CHNS_V &chnsV, CHNS_GFS chnsgfs, int verbosity = 0)
+ {
+ if(verbosity > 1)
+ {
+ printVector(chnsV,verbosity-2,"chnsV");
+ }
+
+ // do refinement
+ if(verbosity > 0) std::cout << "adapt grid and solution... ";
+ my_adapt_grid(grid, gv, chnsgfs, chnsV, 4);
+
+ if(verbosity > 0) std::cout << " ...adapt ended" << std::endl;
+
+ if(verbosity > 1)
+ {
+
+ printVector(chnsV,verbosity-2,"chnsV");
+ }
+ }
+
template<typename Grid, typename CH_V, typename NS_V, typename CH_GFS, typename NS_GFS>
void myRefineGrid(Grid &grid, GV& gv, CH_V &chV, CH_V &chVOld, NS_V &nsV, NS_V &nsVOld, CH_GFS chgfs, NS_GFS nsgfs, int verbosity = 0)
{
@@ -502,9 +601,6 @@ public:
// do refinement
if(verbosity > 0) std::cout << "adapt grid and solution... ";
- //auto chTransfer = transferSolutions(chgfs,4,chV);
- //auto nsTransfer = transferSolutions(nsgfs,4,nsV);
- //Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
my_adapt_grid(grid, gv, chgfs, chV, chVOld, nsgfs, nsV, nsVOld, 4);
if(verbosity > 0) std::cout << " ...adapt ended" << std::endl;
@@ -516,4 +612,24 @@ public:
printVector(nsV,verbosity-2,"nsV");
}
}
+
+ template<typename Grid, typename CHNS_V, typename CHNS_GFS>
+ void myRefineGrid(Grid &grid, GV& gv, CHNS_V &chnsV, CHNS_V &chnsVOld, CHNS_GFS chnsgfs, int verbosity = 0)
+ {
+ if(verbosity > 1)
+ {
+ printVector(chnsV,verbosity-2,"chnsV");
+ }
+
+ // do refinement
+ if(verbosity > 0) std::cout << "adapt grid and solution... ";
+ my_adapt_grid(grid, gv, chnsgfs, chnsV, chnsVOld, 4);
+
+ if(verbosity > 0) std::cout << " ...adapt ended" << std::endl;
+
+ if(verbosity > 1)
+ {
+ printVector(chnsV,verbosity-2,"chnsV");
+ }
+ }
};
diff --git a/dune/phasefield/ffs_chns_r.hh b/dune/phasefield/ffs_chns_r.hh
index f3a386aed7e0a8a0ade8d65f141a45612fea71ae..3e403cdf5236068b03bfb9ab6685615ae6eff71f 100644
--- a/dune/phasefield/ffs_chns_r.hh
+++ b/dune/phasefield/ffs_chns_r.hh
@@ -43,6 +43,7 @@ public:
typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS;
NS ns;
+// typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagB, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS, U_GFS> CH;
typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS, U_GFS> CH;
CH ch;
@@ -79,6 +80,74 @@ public:
}
};
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM, typename U_FEM>
+class GFS_ffs_chns_r_monolithic
+{
+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<VectorBackend, OrderingTagB, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS, U_GFS> CH;
+// typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS, U_GFS> CH;
+ CH ch;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, CH, NS> CHNS;
+ CHNS chns;
+
+ typedef typename CHNS::template ConstraintsContainer<RF>::Type CHNS_CC;
+ CHNS_CC chnsCC;
+
+
+ GFS_ffs_chns_r_monolithic(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),
+ chns(ch,ns),
+ chnsCC()
+ {
+ 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.bdryCHNS,chns,chnsCC, 0);
+ }
+};
+
//_________________________________________________________________________________________________________________________
@@ -162,6 +231,86 @@ public:
}
};
+template<typename GFS, typename CHNS_V>
+class DGF_ffs_chns_r_monolithic
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<1,0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,CHNS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<1,1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,CHNS_V> PDGF;
+ PDGF pDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<0,0> > PhiSubGFS;
+ PhiSubGFS phiSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PhiSubGFS,CHNS_V> PhiDGF;
+ PhiDGF phiDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<0,1> > MuSubGFS;
+ MuSubGFS muSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<MuSubGFS,CHNS_V> MuDGF;
+ MuDGF muDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<0,2> > Phi2SubGFS;
+ Phi2SubGFS phi2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Phi2SubGFS,CHNS_V> Phi2DGF;
+ Phi2DGF phi2Dgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<0,3> > Mu2SubGFS;
+ Mu2SubGFS mu2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Mu2SubGFS,CHNS_V> Mu2DGF;
+ Mu2DGF mu2Dgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CHNS,Dune::TypeTree::TreePath<0,4> > USubGFS;
+ USubGFS uSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<USubGFS,CHNS_V> UDGF;
+ UDGF uDgf;
+
+ //typedef PF_3P_Color_DGF<typename GFS::CH,CH_V,3> ColorDGF;
+ //ColorDGF colorDgf;
+
+
+ DGF_ffs_chns_r_monolithic(const GFS& gfs, const CHNS_V& chnsV) :
+ vSubGfs(gfs.chns),
+ vDgf(vSubGfs,chnsV),
+ pSubGfs(gfs.chns),
+ pDgf(pSubGfs,chnsV),
+
+ phiSubGfs(gfs.chns),
+ phiDgf(phiSubGfs,chnsV),
+ muSubGfs(gfs.chns),
+ muDgf(muSubGfs,chnsV),
+ phi2SubGfs(gfs.chns),
+ phi2Dgf(phi2SubGfs,chnsV),
+ mu2SubGfs(gfs.chns),
+ mu2Dgf(mu2SubGfs,chnsV),
+
+ uSubGfs(gfs.chns),
+ uDgf(uSubGfs,chnsV)
+
+ //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>
@@ -195,6 +344,40 @@ public:
}
};
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu, typename IniPhi2, typename IniMu2, typename IniU>
+class Ini_ffs_chns_r_monolithic
+{
+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;
+ typedef Dune::PDELab::CompositeGridFunction<IniCH,IniNS> IniCHNS;
+ IniCHNS iniCHNS;
+
+ template<typename GV, typename Params>
+ Ini_ffs_chns_r_monolithic(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),
+ iniCHNS(iniCH,iniNS)
+ {
+ }
+};
+
//_________________________________________________________________________________________________________________________
template<typename BdryV1, typename BdryV2, typename BdryP,
@@ -270,6 +453,46 @@ public:
}
};
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2, typename BdryU>
+class Bdry_ffs_chns_r_compositeV_2d_Explicit_Domain_monolithic
+{
+public:
+ BdryV1 bdryV1;
+ BdryV2 bdryV2;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ BdryPhi2 bdryPhi2;
+ BdryMu2 bdryMu2;
+ BdryU bdryU;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu, BdryPhi2, BdryMu2, BdryU> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryCH, BdryNS> BdryCHNS;
+ BdryCHNS bdryCHNS;
+
+ template<typename DomainSize>
+ Bdry_ffs_chns_r_compositeV_2d_Explicit_Domain_monolithic(DomainSize ds) :
+ bdryV1(ds),
+ bdryV2(ds),
+ bdryP(ds),
+ bdryPhi(ds),
+ bdryMu(ds),
+ bdryPhi2(ds),
+ bdryMu2(ds),
+ bdryU(ds),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryPhi2,bdryMu2,bdryU),
+ bdryCHNS(bdryCH,bdryNS)
+ {
+ }
+};
+
template<typename RF>
class DomainSize
{
diff --git a/dune/phasefield/localoperator/chns_monolithic_diff_estimator.hh b/dune/phasefield/localoperator/chns_monolithic_diff_estimator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d03b8f3f2688cd2291caaac0f6b532aabbdf9978
--- /dev/null
+++ b/dune/phasefield/localoperator/chns_monolithic_diff_estimator.hh
@@ -0,0 +1,85 @@
+#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 CHNSGFS>
+class CHNS_MONOLITHIC_DIFF_ESTIMATOR :
+ public Dune::PDELab::LocalOperatorDefaultFlags
+{
+private:
+ typedef typename CHNSGFS::template Child<0>::Type::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
+ CHNS_MONOLITHIC_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& chspace = child(lfsu,_0);
+ const auto& pf1space = child(chspace,_0);
+ const auto& pf2space = child(chspace,_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/dune/phasefield/localoperator/fem_base_chns.hh b/dune/phasefield/localoperator/fem_base_chns.hh
index d01e62e38b75fa21212869684dc39be69d66799a..98de02b985b31252795c331370643fac05ed6bc2 100644
--- a/dune/phasefield/localoperator/fem_base_chns.hh
+++ b/dune/phasefield/localoperator/fem_base_chns.hh
@@ -223,7 +223,7 @@ public:
gradpf2.axpy(lv[pf2space.localIndex(i)],this->basis.gradphi[i][0]);
xi2 += lv[xi2space.localIndex(i)]*this->basis.phi[i];
gradxi2.axpy(lv[xi2space.localIndex(i)],this->basis.gradphi[i][0]);
- }
+ }
}
};
@@ -392,6 +392,106 @@ public:
};
//_________________________________________________________________________________________________________________________
+template<typename RF, const int dim, typename CHNSLFS>
+class CHNSVariables_All
+{
+public:
+ const typename CHNSLFS::template Child<1>::Type& nsspace;
+ const typename CHNSLFS::template Child<1>::Type::template Child<0>::Type& vspace;
+ const typename CHNSLFS::template Child<1>::Type::template Child<0>::Type::template Child<0>::Type& vfemspace;
+ BasisFN<typename CHNSLFS::template Child<1>::Type::template Child<0>::Type::template Child<0>::Type> vbasis;
+ Dune::FieldVector<RF,dim> v;
+ Dune::FieldMatrix<RF,dim,dim> jacv;
+ RF div_v;
+
+ const typename CHNSLFS::template Child<1>::Type::template Child<1>::Type& pspace;
+ BasisFN<typename CHNSLFS::template Child<1>::Type::template Child<1>::Type> pbasis;
+ RF p;
+ Dune::FieldVector<RF,dim> gradp;
+
+ const typename CHNSLFS::template Child<0>::Type& chspace;
+ const typename CHNSLFS::template Child<0>::Type::template Child<0>::Type& pfspace;
+ const typename CHNSLFS::template Child<0>::Type::template Child<1>::Type& xispace;
+ const typename CHNSLFS::template Child<0>::Type::template Child<2>::Type& pf2space;
+ const typename CHNSLFS::template Child<0>::Type::template Child<3>::Type& xi2space;
+ BasisFN<typename CHNSLFS::template Child<0>::Type::template Child<0>::Type> basis;
+ RF pf;
+ RF pf2;
+ RF xi;
+ RF xi2;
+ Dune::FieldVector<RF,dim> gradpf;
+ Dune::FieldVector<RF,dim> gradxi;
+ Dune::FieldVector<RF,dim> gradpf2;
+ Dune::FieldVector<RF,dim> gradxi2;
+
+ const typename CHNSLFS::template Child<0>::Type::template Child<4>::Type& uspace;
+ RF u;
+ Dune::FieldVector<RF,dim> gradu;
+
+
+ template<typename Pos, typename CHCache, typename VCache, typename PCache, typename LV, typename S>
+ CHNSVariables_All(const CHCache& chcache, const VCache& vcache, const PCache& pcache,
+ const Pos& position, const CHNSLFS& chnsLfs, const LV& lv, S& s) :
+ nsspace(child(chnsLfs,Dune::Indices::_1)),
+ vspace(child(nsspace,Dune::Indices::_0)),
+ vfemspace(child(vspace,Dune::Indices::_0)),
+ vbasis(vcache, position, vfemspace, s),
+ v(0.0),
+ jacv(0.0),
+ div_v(0.0),
+ pspace(child(nsspace,Dune::Indices::_1)),
+ pbasis(pcache, position, pspace, s),
+ p(0.0),
+ gradp(0.0),
+ chspace(child(chnsLfs,Dune::Indices::_0)),
+ pfspace(child(chspace,Dune::Indices::_0)),
+ xispace(child(chspace,Dune::Indices::_1)),
+ pf2space(child(chspace,Dune::Indices::_2)),
+ xi2space(child(chspace,Dune::Indices::_3)),
+ basis(chcache, position, pfspace, s),
+ pf(0.0),
+ pf2(0.0),
+ xi(0.0),
+ xi2(0.0),
+ gradpf(0.0),
+ gradxi(0.0),
+ gradpf2(0.0),
+ gradxi2(0.0),
+ uspace(child(chspace,Dune::Indices::_4)),
+ u(0.0),
+ gradu(0.0)
+ {
+ for (size_t i=0; i<pfspace.size(); i++)
+ {
+ pf += lv[pfspace.localIndex(i)]*basis.phi[i];
+ pf2 += lv[pf2space.localIndex(i)]*basis.phi[i];
+ xi += lv[xispace.localIndex(i)]*basis.phi[i];
+ xi2 += lv[xi2space.localIndex(i)]*basis.phi[i];
+ gradpf.axpy(lv[pfspace.localIndex(i)],basis.gradphi[i][0]);
+ gradxi.axpy(lv[xispace.localIndex(i)],basis.gradphi[i][0]);
+ gradpf2.axpy(lv[pf2space.localIndex(i)],basis.gradphi[i][0]);
+ gradxi2.axpy(lv[xi2space.localIndex(i)],basis.gradphi[i][0]);
+ u += lv[uspace.localIndex(i)]*basis.phi[i];
+ gradu.axpy(lv[uspace.localIndex(i)],basis.gradphi[i][0]);
+ }
+ for(size_t k=0; k<dim; k++)
+ {
+ for (size_t i=0; i<vfemspace.size(); i++)
+ {
+ v[k] += lv[child(vspace,k).localIndex(i)]*vbasis.phi[i];
+ jacv[k].axpy(lv[child(vspace,k).localIndex(i)],vbasis.gradphi[i][0]);
+ }
+ div_v += jacv[k][k]; // div is trace of jacobian
+ }
+ for (size_t i=0; i<pspace.size(); i++)
+ {
+ p += lv[pspace.localIndex(i)]*pbasis.phi[i];
+ gradp.axpy(lv[pspace.localIndex(i)],pbasis.gradphi[i][0]);
+ }
+ }
+};
+
+//_________________________________________________________________________________________________________________________
template<typename CHGFS>
using CHLocalBasisCache =
@@ -406,6 +506,19 @@ template<typename NSGFS>
using PLocalBasisCache =
Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<1>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+template<typename CHNSGFS>
+using CH_CHNSLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<0>::Type::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
+
+template<typename CHNSGFS>
+using V_CHNSLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<1>::Type::template Child<0>::Type::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
+template<typename CHNSGFS>
+using P_CHNSLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename CHNSGFS::template Child<1>::Type::template Child<1>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
+
// template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
// class FEM_CHNS_Base :
// public Dune::PDELab::FullVolumePattern,
diff --git a/dune/phasefield/porenetwork/1p_chns.hh b/dune/phasefield/porenetwork/1p_chns.hh
new file mode 100644
index 0000000000000000000000000000000000000000..843dc73aa4bf679e537aff9af66cb765d9b02ac1
--- /dev/null
+++ b/dune/phasefield/porenetwork/1p_chns.hh
@@ -0,0 +1,160 @@
+#pragma once
+
+
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/subspace.hh>
+#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename RF, typename ES, typename V_FEM, typename PHI_FEM>
+class GFS_1p_chns
+{
+public:
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
+ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
+ typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
+ typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
+
+ ConstraintsAssembler conphi, conmu;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, V_FEM, ConstraintsAssembler, VectorBackend> VX_GFS;
+ VX_GFS vxGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> PHI_GFS;
+ PHI_GFS phiGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> MU_GFS;
+ MU_GFS muGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, VX_GFS> NS;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH;
+ CH ch;
+
+ typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
+ NS_CC nsCC;
+
+ typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
+ CH_CC chCC;
+
+ GFS_1p_chns(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
+ conphi(),conmu(),
+ vxGfs(es,vFem),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ ns(vxGfs),
+ ch(phiGfs, muGfs),
+ nsCC(),
+ chCC()
+ {
+ vxGfs.name("velocity");
+ phiGfs.name("phasefield");
+ muGfs.name("potential");
+ }
+
+ template<typename Bdry>
+ void updateConstraintsContainer(const Bdry& bdry)
+ {
+ Dune::PDELab::constraints(bdry.bdryNS,ns,nsCC, 0);
+ Dune::PDELab::constraints(bdry.bdryCH,ch,chCC, 0);
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename GFS, typename NS_V, typename CH_V>
+class DGF_1p_chns
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VXSubGFS;
+ VXSubGFS vxSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<VXSubGFS,NS_V> VXDGF;
+ VXDGF vxDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<0> > PhiSubGFS;
+ PhiSubGFS phiSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PhiSubGFS,CH_V> PhiDGF;
+ PhiDGF phiDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<1> > MuSubGFS;
+ MuSubGFS muSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<MuSubGFS,CH_V> MuDGF;
+ MuDGF muDgf;
+
+
+ DGF_1p_chns(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
+ vxSubGfs(gfs.ns),
+ vxDgf(vxSubGfs,nsV),
+ phiSubGfs(gfs.ch),
+ phiDgf(phiSubGfs,chV),
+ muSubGfs(gfs.ch),
+ muDgf(muSubGfs,chV)
+ {
+ }
+
+ template<typename VTKWRITER>
+ void addToVTKWriter(VTKWRITER& vtkwriter)
+ {
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vxDgf,"vx"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename IniVX, typename IniPhi, typename IniMu>
+class Ini_1p_chns
+{
+public:
+ IniVX iniVx;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ typedef Dune::PDELab::CompositeGridFunction<IniVX> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ template<typename GV, typename Params>
+ Ini_1p_chns(const GV& gv, const Params& params) :
+ iniVx(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniNS(iniVx),
+ iniCH(iniPhi,iniMu)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename BdryVX, typename BdryPhi, typename BdryMu>
+class Bdry_1p_chns
+{
+public:
+ BdryVX bdryVx;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ Bdry_ff_chns() :
+ bdryVx(),
+ bdryPhi(),
+ bdryMu(),
+ bdryNS(bdryVx),
+ bdryCH(bdryPhi,bdryMu)
+ {
+ }
+};
diff --git a/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh b/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48158406b535fe449ffa2beab933767156d5afe8
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
@@ -0,0 +1,88 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+//Parameter class (for functions), FEM space and RangeField
+template<typename Param, typename CHGFS, typename CHContainer>
+class FEM_1P_CahnHilliard_R :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_1P_CahnHilliard_R<Param, CHGFS, CHContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_1P_CahnHilliard_R<Param, CHGFS, CHContainer> >
+{
+private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chOldLocal;
+
+ Param& param; // parameter functions
+
+public:
+ enum {doPatternVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_1P_CahnHilliard_R(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) :
+ chOldLocal(chGfs_,chVOld_),
+ param(param_)
+ {
+ }
+
+
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+ {
+ //Quadrature Rule
+ const int dim = EG::Entity::dimension;
+ const int order = 4; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chOldLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(), chOldLocal.lv, S);
+
+
+
+ for (size_t i=0; i<ch.pfspace.size(); i++)
+ {
+
+ //RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
+ RF Rspeed = param.dw.eval_q(ch.pfOld);
+ RF R1 = -Rspeed/param.eps*param.reaction.eval(1)*ch.basis.phi[i];
+
+ RF J1 = param.mobility.eval(ch.pf) * (ch.gradxi*ch.basis.gradphi[i][0]);
+
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ +J1-R1
+ ));
+
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.dw.eval_d(ch.pf) *ch.basis.phi[i]
+ - param.eps * (ch.gradpf*ch.basis.gradphi[i][0])
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d6f9656ff41726467b8084de1b1cc74040a0a6c9
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
@@ -0,0 +1,92 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+#include"fem_base_porenetwork.hh"
+
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+ class FEM_1P_NavierStokes :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+ {
+ private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+ Micro_VXLocalBasisCache<NSGFS> vCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chLocal;
+
+ Param& param; // parameter functions
+
+ //FOR DETECTING ENTITY changes
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
+ mutable ENT OldEnt;
+
+ public:
+
+
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ TMIXED_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
+ NSGFS& nsGfs_) :
+ chLocal(chGfs_,chV_),
+ param(param_)
+ {
+ }
+
+
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+ {
+// std::cout << "hi_alpha" << std::endl;
+ //Quadrature Rule
+ const int dim = EG::Entity::dimension;
+ const int order = 5; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_Threephase_Micro_Extension<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+ NSMicroVariables<RF,dim,LFSV> ns(vCache, ip.position(), lfsv, x.base(), S);
+ RF pf3 = 1-ch.pf - ch.pf2;
+ RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+
+ for (size_t i=0; i<ns.vxspace.size(); i++)
+ {
+ r.accumulate(ns.vxspace, i, factor*(
+ param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]
+ + param.phys.eval_mu(pf_f) * (ns.gradvx * ns.vbasis.gradphi[i][0])
+ - pf_f * ns.vbasis.phi[i]
+ ));
+ } // for i
+ } // for ip
+// std::cout << "bye_alpha" << std::endl;
+ } //alpha volume
+
+ };
diff --git a/dune/phasefield/porenetwork/fem_base_porenetwork.hh b/dune/phasefield/porenetwork/fem_base_porenetwork.hh
new file mode 100644
index 0000000000000000000000000000000000000000..601757c01d5dad6338178ae3e25fd0d5c0a0810f
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_base_porenetwork.hh
@@ -0,0 +1,34 @@
+#pragma once
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+template<typename RF, const int dim, typename NSLFS>
+class NSPorenetworkFlow
+{
+public:
+ const typename NSLFS::template Child<0>::Type& vxspace;
+ BasisFN<typename NSLFS::template Child<0>::Type> vbasis;
+
+ RF vx;
+ Dune::FieldVector<RF,dim> gradvx;
+
+ template<typename Pos, typename vCache, typename LV, typename S>
+ NSPorenetworkFlow(const vCache& vcache, const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
+ vxspace(child(nsLfs,Dune::Indices::_0)),
+ vbasis(vcache, position, vxspace, s),
+ vx(0.0),
+ gradvx(0.0)
+ {
+ for (size_t i=0; i<this->vxspace.size(); i++)
+ {
+ vx += lv[vxspace.localIndex(i)]*vbasis.phi[i];
+ gradvx.axpy(lv[vxspace.localIndex(i)],vbasis.gradphi[i][0]);
+ }
+ }
+};
+
+
+template<typename NSGFS>
+using Micro_VXLocalBasisCache =
+ Dune::PDELab::LocalBasisCache<typename NSGFS::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..71f32f12c91fd5091c3a063a4b55fb963847c2af
--- /dev/null
+++ b/dune/phasefield/porenetwork/pn_1pflow.hh
@@ -0,0 +1,567 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+#include<vector>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include <dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh>
+#include "tmixed_chns_navierstokes.hh"
+#include "tmixed_upscaledp.hh"
+#include "tmixed_upscaledc.hh"
+#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include "tmixed_parameters.hh"
+#include "tmixed_boundary.hh"
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/porenetwork/1p_chns.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+template< typename Grid>
+std::shared_ptr<Grid> initPnGrid(std::string filename)
+{
+ Dune::GridFactory<Grid> factory;
+ Dune::GmshReader<Grid>::read(factory,filename,true,true);
+ return std::shared_ptr<Grid>(factory.createGrid());
+}
+
+
+template<typename Parameters>
+class Pn_1PFlow{
+public:
+
+ typedef double RF;
+ const int dim = 2;
+
+ typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+ typedef Grid::ctype DF;
+ typedef Grid::LeafGridView GV;
+ typedef Dune::PDELab::AllEntitySet<GV> ES;
+ const int k = 2;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+ const int l = 1;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+ typedef GFS_1p_chns<RF, ES, V_FEM, PHI_FEM> GFS;
+
+ typedef typename GFS::NS NS_GFS;
+ typedef typename GFS::CH CH_GFS;
+ typedef Dune::PDELab::Backend::Vector<CH_GFS, RF> CH_V;
+ typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
+
+ typedef Ini_1p_chns< ZeroInitial<GV,RF,Parameters>, //V
+ PhiInitial_Thinstrip_Layers<GV,RF,Parameters,0>, //TODO
+ ZeroInitial<GV,RF,Parameters>> Initial;
+
+ typedef Bdry_1p_chns<Dune::PDELab::DirichletConstraintsParameters, //V //TODO
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
+
+ typedef TMIXED_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+ typedef FEM_1P_CahnHilliard_R<Parameters, CH_GFS, CH_V> CH_LOP;
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+
+ typedef DGF_1p_chns<GFS, NS_V, CH_V> DGF;
+
+ typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+ typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS;
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS;
+
+ typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> NS_Newton;
+ typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> CH_Newton;
+
+
+ std::shared_ptr<Grid> gridp;
+ GV gv;
+ ES es;
+ V_FEM vFem;
+ P_Fem pFem;
+ PHI_FEM phiFem;
+
+ GFS gfs;
+
+ CH_V chV, chVOld;
+ NS_V nsV;
+ Initial initial;
+ Bdry bdry;
+
+
+ RF dxp = 1;
+ RF c = 1;
+ MBE mbe;
+
+ RF kf = 0;
+ RF phic = 0;
+ RF phicOld = 0;
+ RF phiSolid = 0;
+ RF phiSolidOld = 0;
+ RF kc = 0;
+
+ std::unique_ptr<NS_LOP> nsLop;
+ std::unique_ptr<CH_LOP> chLop;
+ std::unique_ptr<DGF> dgf;
+ std::unique_ptr<NS_GO> nsGo;
+ std::unique_ptr<CH_GO> chGo;
+ std::unique_ptr<NS_LS> nsLs;
+ std::unique_ptr<CH_LS> chLs;
+ std::unique_ptr<NS_Newton> nsNewton;
+ std::unique_ptr<CH_Newton> chNewton;
+
+ Pn_1PFlow(std::string filename, Parameters& param) :
+ gridp(initPnGrid<Grid>(filename)),
+ gv(gridp->leafGridView()),
+ es(gv),
+ vFem(es), phiFem(es),
+ gfs(es,vFem,phiFem),
+ chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch),
+ initial(gv,param), mbe(10),
+ {
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+
+ chVOld = chV;
+
+ //Local Operator
+ nsLop = std::make_unique<NS_LOP>(param, gfs.ch, chV, gfs.ns);
+ chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld);
+
+ dgf = std::make_unique<DGF>(gfs,nsV,chV);
+ }
+
+ void initializeGridOperators(Parameters& param)
+ {
+ //Grid Operator
+ nsGo = std::make_unique<NS_GO>(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nsLop,mbe);
+ chGo = std::make_unique<CH_GO>(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,*chLop,mbe);
+ // Linear solver
+ chLs = std::make_unique<CH_LS>(*chGo,500,100,3,1);
+ nsLs = std::make_unique<NS_LS>(*nsGo,1000,200,3,1);
+ // Nonlinear solver
+
+ nsNewton = std::make_unique<NS_Newton>(*nsGo,nsV,*nsLs);
+ nsNewton->setParameters(param.nsNewtonTree);
+
+ chNewton = std::make_unique<CH_Newton>(*chGo,chV,*chLs);
+ chNewton->setParameters(param.chNewtonTree);
+
+ printVector(chV,10,"chV");
+ printVector(nsV,10,"nsV");
+
+ updatePhiC();
+ updatePhiSolid();
+ }
+
+ void applyNs()
+ {
+ int verb = 0;
+ if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
+ nsNewton->apply();
+ nsVOld = nsV;
+ updateKf();
+ updateKc();
+ }
+
+ void applyCh()
+ {
+ phicOld = phic;
+ phiSolidOld = phiSolid;
+ int verb=0;
+ if (verb > 1) std::cout << "= Begin newtonapply CH:" << " with dxp= " << dxp << " dxpUpwind= " << dxpUpwind << " c= " << c << std::endl;
+ chNewton->apply();
+ chVOld = chV;
+ updatePhiC();
+ updatePhiSolid();
+ }
+
+ void updateKf()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> w;
+ dgfMicro->vxDgf.evaluate(element, xlocal, w);
+ Dune::FieldVector<double,1> phi1;
+ dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ Dune::FieldVector<double,1> phi2;
+ dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+ return Dune::FieldVector<double,1>((phi1+phi2)*w);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ kf = integral[0];
+ }
+
+ void updatePhiC()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> phi1;
+ dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ return Dune::FieldVector<double,1>(phi1);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ phic = integral[0];
+ }
+
+ void updatePhiSolid()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> phi1;
+ dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ Dune::FieldVector<double,1> phi2;
+ dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+ return Dune::FieldVector<double,1>(1-phi1-phi2);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ phiSolid = integral[0];
+ }
+
+ void updateKc()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> w;
+ dgfMicro->vxDgf.evaluate(element, xlocal, w);
+ Dune::FieldVector<double,1> phi1;
+ dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ return Dune::FieldVector<double,1>(phi1*w);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ kc = integral[0];
+ }
+};
+
+template<typename RF, typename YSimS>
+class P_ParamAdaptor
+{
+public:
+ const int xs;
+ YSimS& ySimS;
+ P_ParamAdaptor(YSimS& ySimS_) : ySimS(ySimS_), xs(ySimS_.size())
+ {
+ }
+
+ RF evalKf(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].kf;
+ }
+
+ RF evalKc(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].kc;
+ }
+
+ RF evalDxP(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].dxp;
+ }
+
+
+ RF evalPhiC(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].phic;
+ }
+
+ RF evalPhiCOld(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].phicOld;
+ }
+
+ RF evalPhiSolid(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].phiSolid;
+ }
+
+ RF evalPhiSolidOld(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].phiSolidOld;
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template< typename XGV, typename PGrad_DGF, typename SimS>
+void evaluateDxP(XGV& xgv, PGrad_DGF& pGradDgf, SimS& simS)
+{
+ Dune::FieldVector<double,1> pgrad(0);
+ Dune::FieldVector<double,1> pgradUpwind(0);
+ int xs = simS.size();
+ auto it = elements(xgv).begin();
+ for(int i=0; i<xs; i++)
+ {
+ pGradDgf.evaluate(*it,it->geometry().local(it->geometry().center()), pgrad );
+ pgrad[0] = pgrad[0] - 1; //TODO make more general with pressure drop/distance
+ simS[i].dxp = pgrad[0];
+ simS[i].dxpUpwind = pgradUpwind[0];
+ it++;
+ pgradUpwind = pgrad;
+ }
+ //Periodic:
+ simS[0].dxpUpwind = simS[xs-1].dxp;
+}
+
+template< typename XGV, typename C_DGF, typename SimS>
+void evaluateC(XGV& xgv, C_DGF& cDgf, SimS& simS)
+{
+ Dune::FieldVector<double,1> cTmp(0);
+ int xs = simS.size();
+ auto it = elements(xgv).begin();
+ for(int i=0; i<xs; i++)
+ {
+ cDgf.evaluate(*it,it->geometry().local(it->geometry().center()), cTmp );
+ simS[i].c = cTmp[0];
+ it++;
+ }
+}
+
+//_________________________________________________________________________________________________________________________
+
+template<typename XGV, typename YGV, typename X_ES, typename Y_ES, typename TGV, typename T_FEM, typename V_FEM, typename P_FEM, typename PHI_FEM, typename C_FEM>
+void driver_mixed(XGV& xgv, X_ES& x_es, const P_FEM& pFem, const C_FEM& cFem, const int xsize,
+ YGV& ygv, Y_ES& y_es, TGV& tgv, const T_FEM& tFem, const V_FEM& vFem, const PHI_FEM& phiFem,
+ Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+ const int dim = 1;
+ typedef double RF;
+
+
+ typedef Params_tmixed< RF,
+ //TripleWell<RF,DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> > >,
+ TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >,
+ Mobility_Const<RF>,
+ Reaction_Thinstrip<RF>,
+ VelDissipation_Quadratic<RF>
+ //VelDissipation_Quadratic_Shifted<RF>
+ > Parameters;
+ Parameters param(ptree);
+ TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make micro problems
+
+
+ typedef GFS_micro<RF, Y_ES, V_FEM, PHI_FEM> GFS;
+ GFS gfs(y_es,vFem,phiFem);
+
+ typedef YSlice<YGV, GFS, Parameters> YSim;
+ typedef std::vector<YSim> YSimS;
+ YSimS ySimS;
+ ySimS.reserve(xsize);
+ for(int i=0; i<xsize; i++)
+ {
+ ySimS.emplace_back(ygv, gfs, param, ((RF)i)/((RF)xsize), ((RF)1)/((RF)xsize) );
+ }
+
+
+ //_________________________________________________________________________________________________________________________
+ // Make macro problems
+ typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::NoConstraints NCon;
+ ConstraintsAssembler conP();
+ typedef Dune::PDELab::GridFunctionSpace<X_ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs(x_es,pFem);
+
+ auto blambda = [](const auto& x){return true;};
+ auto b = Dune::PDELab::
+ makeBoundaryConditionFromCallable(xgv,blambda);
+ //typedef Dune::PDELab::DirichletConstraintsParameters BdryP;
+ //BdryP pBdry();
+ typedef typename P_GFS::template ConstraintsContainer<RF>::Type P_CC;
+ P_CC pCC;
+ //pCC.clear();
+ Dune::PDELab::constraints(b,pGfs,pCC);
+
+ typedef Dune::PDELab::Backend::Vector<P_GFS, RF> P_V;
+ P_V pV(pGfs);
+
+ auto p0lambda = [&](const auto& xg)
+ {return (1-xg[0])*(1-xg[0]);};
+ auto p0 = Dune::PDELab::
+ makeGridFunctionFromCallable(xgv,p0lambda);
+ Dune::PDELab::interpolate(p0,pGfs,pV);
+
+
+
+
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+ MBE mbe(10);
+ typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
+ P_PAdaptor pPAdaptor(ySimS);
+ typedef TMIXED_UpscaledP<P_PAdaptor,P_FEM,RF> P_LOP;
+ P_LOP pLop(pPAdaptor);
+
+ typedef Dune::PDELab::GridOperator<P_GFS, P_GFS, P_LOP, MBE,
+ RF, RF, RF, P_CC, P_CC> P_GO;
+ P_GO pGo(pGfs, pCC, pGfs, pCC, pLop, mbe);
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<P_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> P_LS;
+ P_LS pLs(pGo,500,100,3,1);
+ typedef Dune::PDELab::Newton<P_GO,P_LS,P_V> P_Newton;
+ P_Newton pNewton(pGo,pV,pLs);
+ pNewton.setReassembleThreshold(0.0); // always reassemble J
+ pNewton.setVerbosityLevel(1); // be verbose
+ pNewton.setReduction(1e-10); // total reduction
+ pNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
+ pNewton.setMaxIterations(25); // limit number of its
+ pNewton.setLineSearchMaxIterations(10); // limit line search
+ pNewton.setAbsoluteLimit(1e-10);
+
+
+
+
+
+
+ typedef Dune::PDELab::GridFunctionSpace<X_ES, C_FEM,Dune::PDELab::NoConstraints,VectorBackend> C_GFS;
+ C_GFS cGfs(x_es,cFem);
+ typedef Dune::PDELab::Backend::Vector<C_GFS, RF> C_V;
+ C_V cV(cGfs);
+ C_V cVOld(cGfs);
+ auto c0lambda = [&](const auto& xg)
+ //{return 0.5+0.4*std::sin(4*3.14152*xg[0]);};
+ {return 0.5;};
+ auto c0 = Dune::PDELab::
+ makeGridFunctionFromCallable(xgv,c0lambda);
+ Dune::PDELab::interpolate(c0,cGfs,cV);
+ cVOld = cV;
+ //typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
+ //P_PAdaptor pPAdaptor(ySimS);
+ typedef TMIXED_UpscaledC<Parameters,P_PAdaptor,C_GFS,C_V,RF> C_LOP;
+ C_LOP cLop(param,pPAdaptor,cGfs,cVOld);
+
+ typedef Dune::PDELab::EmptyTransformation NoCC;
+ //typedef typename C_GFS::template ConstraintsContainer<RF>::Type C_CC;
+ //C_CC cCC;
+
+ typedef Dune::PDELab::GridOperator<C_GFS, C_GFS, C_LOP, MBE,
+ RF, RF, RF, NoCC, NoCC> C_GO;
+ MBE mbe2(10);
+ C_GO cGo(cGfs, cGfs, cLop, mbe);
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<C_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> C_LS;
+ C_LS cLs(cGo,500,100,3,1);
+ //typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<C_GO> C_LS;
+ //C_LS cLs(100,2);
+
+ typedef Dune::PDELab::Newton<C_GO,C_LS,C_V> C_Newton;
+ C_Newton cNewton(cGo,cV,cLs);
+ cNewton.setReassembleThreshold(0.0); // always reassemble J
+ cNewton.setVerbosityLevel(1); // be verbose
+ cNewton.setReduction(1e-10); // total reduction
+ cNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
+ cNewton.setMaxIterations(25); // limit number of its
+ cNewton.setLineSearchMaxIterations(10); // limit line search
+ cNewton.setAbsoluteLimit(1e-10);
+
+
+
+
+
+
+ //_________________________________________________________________________________________________________________________
+ // prepare VTK writer
+ std::string filename=ptree.get("output.filename","output");
+ std::string x_filename = "xoutput";
+ 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;
+ }
+ if( stat( x_filename.c_str(), &st ) != 0 )
+ {
+ int stat = 0;
+ stat = mkdir(x_filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << x_filename << std::endl;
+ }
+ auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<TGV> >(tgv,0);
+ Dune::VTKSequenceWriter<TGV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+
+ DGF_macro<RF, YGV, TGV, YSimS> dgfMacro(ygv, tgv,ySimS);
+ dgfMacro.addToVTKWriter(vtkwriter);
+ vtkwriter.write(param.time.t,Dune::VTK::base64);
+
+ auto x_stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<XGV> >(xgv,0);
+ Dune::VTKSequenceWriter<XGV> x_vtkwriter(x_stationaryvtkwriter,x_filename.c_str(),x_filename.c_str(),"");
+
+ typedef Dune::PDELab::DiscreteGridFunction<P_GFS,P_V> P_DGF;
+ P_DGF pDgf(pGfs,pV);
+ typedef Dune::PDELab::DiscreteGridFunctionGradient<P_GFS,P_V> PGrad_DGF;
+ PGrad_DGF pGradDgf(pGfs,pV);
+ typedef Dune::PDELab::DiscreteGridFunction<C_GFS,C_V> C_DGF;
+ C_DGF cDgf(cGfs,cV);
+
+ typedef Dune::PDELab::VTKGridFunctionAdapter<P_DGF> P_VTKF;
+ typedef Dune::PDELab::VTKGridFunctionAdapter<C_DGF> C_VTKF;
+ x_vtkwriter.addVertexData(std::shared_ptr<P_VTKF>(
+ new P_VTKF(pDgf,"p")));
+ x_vtkwriter.addVertexData(std::shared_ptr<C_VTKF>(
+ new C_VTKF(cDgf,"c")));
+ x_vtkwriter.write(param.time.t,Dune::VTK::base64);
+
+//_________________________________________________________________________________________________________________________
+while (!timestepManager.isFinished())
+{
+ for(int i=1; i<xsize; i++)
+ {
+ ySimS[i].chVUpwind = ySimS[i-1].chV;
+ }
+ ySimS[0].chVUpwind = ySimS[xsize-1].chV; //Periodic
+
+ for(int i=0; i<xsize; i++)
+ {
+ ySimS[i].applyNs();
+ }
+
+ for(int i=1; i<xsize; i++)
+ {
+ ySimS[i].nsVUpwind = ySimS[i-1].nsV;
+ }
+ ySimS[0].nsVUpwind = ySimS[xsize-1].nsV; //Periodic
+
+
+ //For Debug
+ //timestepManager.applyTimestep(vtkwriter);
+ //x_vtkwriter.write(param.time.t);
+
+
+ std::cout << "Begin newtonapply Pressure:" << std::endl;
+ pNewton.apply();
+ evaluateDxP(xgv, pGradDgf,ySimS);
+
+ for(int i=0; i<xsize; i++)
+ {
+ ySimS[i].applyCh();
+ }
+ //while (!timestepManager.isFinished())
+ //{
+ std::cout << "Begin newtonapply Concentration:" << std::endl;
+ cNewton.apply();
+ cVOld =cV;
+ evaluateC(xgv,cDgf,ySimS);
+
+
+ timestepManager.applyTimestep(vtkwriter, x_vtkwriter);
+ }
+
+}
diff --git a/src/paper_thinstrip_full/3p_driver_monolithic.hh b/src/paper_thinstrip_full/3p_driver_monolithic.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4f1b9de686f8703addbbfe2897ba83e445ecc21c
--- /dev/null
+++ b/src/paper_thinstrip_full/3p_driver_monolithic.hh
@@ -0,0 +1,206 @@
+#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 "thinstrip_full_chns_monolithic.hh"
+#include <dune/phasefield/localoperator/chns_monolithic_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include "thinstrip_parameters.hh"
+#include "thinstrip_boundary.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 GV, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_3p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+ const int dim = GV::dimension;
+ typedef double RF;
+
+
+ typedef Params_thinstrip< RF,
+ //TripleWell<RF,DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> > >,
+ TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >,
+ Mobility_Const<RF>,
+ Reaction_Thinstrip<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_monolithic<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+ GFS gfs(es,vFem,pFem,phiFem,phiFem);
+ typedef typename GFS::CHNS CHNS_GFS;
+
+
+ //_________________________________________________________________________________________________________________________
+ //Constraints
+ typedef DomainSize<RF> DS;
+ DS dS(0,1,0,ptree.get("domain.len",1),1e-8);
+ typedef Bdry_ffs_chns_r_compositeV_2d_Explicit_Domain_monolithic<
+ BottomBoundary<DS>, // v_x
+ AllBoundary<DS>, // v_y
+ NoBoundary<DS>, // p
+ NoBoundary<DS>, // phi1
+ NoBoundary<DS>, // mu1
+ NoBoundary<DS>, // phi2
+ NoBoundary<DS>, // mu2
+ NoBoundary<DS>> // u
+ Bdry;
+ Bdry bdry(dS);
+ gfs.updateConstraintsContainer(bdry);
+
+ //_________________________________________________________________________________________________________________________
+ //Coefficient vectors
+
+ using CHNS_V = Dune::PDELab::Backend::Vector<CHNS_GFS, RF>;
+ CHNS_V chnsV(gfs.chns);
+
+ //_________________________________________________________________________________________________________________________
+ //Initial Values
+
+ typedef Ini_ffs_chns_r_monolithic< VelocityInitial_InOutFlow<GV,RF,Parameters,dim>,
+ PInitial_InOutFlow<GV,RF,Parameters>,
+ PhiInitial_Thinstrip_Layers<GV,RF,Parameters,0>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_Thinstrip_Layers<GV,RF,Parameters,1>,
+ ZeroInitial<GV,RF,Parameters>,
+ UInitial_Inflow<GV,RF,Parameters> > Initial;
+ Initial initial(gv,param);
+ Dune::PDELab::interpolate(initial.iniCHNS, gfs.chns, chnsV);
+
+ //_________________________________________________________________________________________________________________________
+// Refine Intial Values
+ for(int i=0; i<ptree.get("RefinementInitial.maxRefineLevel",2); i++)
+ {
+ // mark elements for refinement
+
+
+ typedef RefinementIndicator<GV,RF,CHNS_MONOLITHIC_DIFF_ESTIMATOR<CHNS_GFS>,CHNS_GFS> RIND0;
+ RIND0 rInd0(gv,gfs.chns,ptree.sub("RefinementInitial"));
+ rInd0.calculateIndicator(gv,chnsV,false);
+ rInd0.markGrid(grid,2);
+ rInd0.myRefineGrid(grid,gv,chnsV,gfs.chns,1);
+ Dune::PDELab::interpolate(initial.iniCHNS , gfs.chns, chnsV);
+ }
+
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry_monolithic(gfs,initial.iniCHNS,chnsV);
+
+ CHNS_V chnsVOld = chnsV;
+
+ //_________________________________________________________________________________________________________________________
+ // LocalOperators
+
+
+
+ typedef THINSTRIP_FULL_CHNS_Monolithic<Parameters, CHNS_GFS, CHNS_V> CHNS_LOP;
+ CHNS_LOP chnsLop(param, gfs.chns, chnsVOld);
+
+ 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_monolithic<GFS, CHNS_V> DGF;
+ DGF dgf(gfs, chnsV);
+
+ //_________________________________________________________________________________________________________________________
+ // 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(chnsV,20,"chnsV");
+
+
+ //____________________________________________________________________________________________________________________
+ //int i_it=0;
+ chnsVOld = chnsV;
+ while (!timestepManager.isFinished())
+ {
+
+ // Grid Operators
+ typedef Dune::PDELab::GridOperator<CHNS_GFS,CHNS_GFS,CHNS_LOP,MBE,RF,RF,RF,typename GFS::CHNS_CC,typename GFS::CHNS_CC> CHNS_GO;
+ CHNS_GO chnsGo(gfs.chns,gfs.chnsCC,gfs.chns,gfs.chnsCC,chnsLop,mbe);
+
+ // Linear solver
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+ <CHNS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CHNS_LS; CHNS_LS ls_chns(chnsGo,1500,200,3,1);
+ //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 CHNS_LS; CHNS_LS ls_chns(200,5000,1);
+
+ // Nonlinear solver
+ typedef Dune::PDELab::Newton<CHNS_GO,CHNS_LS,CHNS_V> PDESolver1;
+ PDESolver1 chns_newton(chnsGo,chnsV,ls_chns);
+ chns_newton.setParameters(ptree.sub("NS_Newton"));
+
+ //_________________________________________________________________________________________________________________________
+ // do time step
+
+
+ interpolateToBdry_monolithic(gfs,initial.iniCHNS, chnsV);
+
+ int result = timestepManager.doOnlyNSTimestep(chns_newton,vtkwriter);
+ if(result <0)
+ {
+ return;
+ }
+ else if(result == 0)
+ {
+ chnsV = chnsVOld;
+ continue;
+ }
+
+ chnsVOld = chnsV;
+//TODO?
+//child(initial.iniCH,Dune::Indices::_0).ApplyReactionDt(param);
+ //_________________________________________________________________________________________________________________________
+
+ // mark elements for refinement
+
+ typedef RefinementIndicator<GV,RF,CHNS_MONOLITHIC_DIFF_ESTIMATOR<CHNS_GFS>,CHNS_GFS> RIND;
+ RIND rInd(gv,gfs.chns,ptree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chnsV,false);
+ rInd.markGrid(grid,2);
+ rInd.myRefineGrid(grid,gv,chnsV,chnsVOld,gfs.chns,1);
+
+ vtkwriter.write(param.time.t+1e-6,Dune::VTK::ascii);
+
+ // recompute constraints
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry_monolithic(gfs,initial.iniCHNS,chnsV);
+
+ //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+ }
+}
diff --git a/src/paper_thinstrip_full/chns.ini b/src/paper_thinstrip_full/chns.ini
index c0f2d28d0006d263c465cc477f99d23e68fb4285..1171fc4692773fb7051eb247ad091fc56f24bb30 100644
--- a/src/paper_thinstrip_full/chns.ini
+++ b/src/paper_thinstrip_full/chns.ini
@@ -22,9 +22,9 @@ YScaling = 0.5
filename = output4
[Phasefield]
-eps=0.03
-delta=0.03
-DoubleWellDelta=0.05
+eps=0.03 #0.03
+delta=0.03 #0.03
+DoubleWellDelta=0.05 #0.05
Mobility=0.01
Sigma1=1
Sigma2=1
@@ -32,23 +32,29 @@ Sigma3=1
[Time]
tMax = 300
-dt = 0.015
+dt = 0.0015
dtMax = 0.15
-dtMin = 0.0015
+dtMin = 0.00015
#saveintervall = 0.5
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
[Initial]
# FlowVyy= 10# pConst/(Domain_width * mu)
-pConst=0.1
-uConst=0.5
-uInflow=0.3
+pConst=0.1 #0.1
+uConst=0.5 #0.5
+uInflow=0.3 #0.3
+
+[RefinementInitial]
+etaRefine = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 4
+minRefineLevel = -1
[Refinement]
etaRefine = 0.2
etaCoarsen = 0.1
-maxRefineLevel = 6
+maxRefineLevel = 4
minRefineLevel = -1
[NS_Newton]
@@ -57,10 +63,12 @@ LineSearchMaxIterations = 30
MaxIterations = 30
AbsoluteLimit = 1e-9
Reduction = 1e-7
-VerbosityLevel = 1
+VerbosityLevel = 2
[Phasefield_Newton]
ReassembleThreshold = 0.0
LineSearchMaxIterations = 30
-MaxIterations = 30
-VerbosityLevel = 1
+MaxIterations = 100
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 2
diff --git a/src/paper_thinstrip_full/thinstrip_boundary.hh b/src/paper_thinstrip_full/thinstrip_boundary.hh
index e241f05b0c31516fa8f5f2cc443f0361bdbf4f34..4ae02510a0896c258289578d7040832e7459a91d 100644
--- a/src/paper_thinstrip_full/thinstrip_boundary.hh
+++ b/src/paper_thinstrip_full/thinstrip_boundary.hh
@@ -192,13 +192,15 @@ public:
{
eps = params.eps;
dSolid = params.initial.ptree.get("dSolid",0.3);
- d1 = params.initial.ptree.get("d1",0.3);
+ d1 = params.initial.ptree.get("d1",0.4);
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
- RF rf = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid));
- RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid - d1));
+ RF x_nondim = x[1]/pref.YScaling;
+ RF d_sine = 0.15*sin(x[0]*2*M_PI);
+ RF rf = pref.tw.dw.eval_shape(1/eps*(x_nondim - dSolid));
+ RF rf2 = pref.tw.dw.eval_shape(1/eps*(x_nondim - dSolid - d1-d_sine));
if(Phase == 1)
{
y = rf2;
diff --git a/src/paper_thinstrip_full/thinstrip_equation.cc b/src/paper_thinstrip_full/thinstrip_equation.cc
index 60d1d9bdbe8f17da7bb0d710baec90a9d99fcebd..0b751ff434ea18a722b7ef59d3108192563d7aff 100644
--- a/src/paper_thinstrip_full/thinstrip_equation.cc
+++ b/src/paper_thinstrip_full/thinstrip_equation.cc
@@ -19,7 +19,8 @@
// pdelab includes
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include <dune/xt/grid/view/periodic.hh>
-#include "3p_driver.hh"
+#include "3p_driver_monolithic.hh"
+//#include "3p_driver.hh"
//===============================================================
// Main program with grid setup
diff --git a/src/paper_thinstrip_full/thinstrip_full_chns_monolithic.hh b/src/paper_thinstrip_full/thinstrip_full_chns_monolithic.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0373565a4a0a060982ecadb68d6e9a64fa57ee60
--- /dev/null
+++ b/src/paper_thinstrip_full/thinstrip_full_chns_monolithic.hh
@@ -0,0 +1,211 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+//Parameter class (for functions), FEM space and RangeField
+template<typename Param, typename CHNSGFS, typename CHNSContainer>
+class THINSTRIP_FULL_CHNS_Monolithic :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<THINSTRIP_FULL_CHNS_Monolithic<Param, CHNSGFS, CHNSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<THINSTRIP_FULL_CHNS_Monolithic<Param, CHNSGFS, CHNSContainer> >
+{
+private:
+ typedef typename CHNSContainer::field_type RF;
+ typedef typename CHNSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CH_CHNSLocalBasisCache<CHNSGFS> chCache;
+ V_CHNSLocalBasisCache<CHNSGFS> vCache;
+ P_CHNSLocalBasisCache<CHNSGFS> pCache;
+
+ typedef GFS_LocalViewer<CHNSGFS, CHNSContainer, ENTITY, RF> CHNSLView;
+ mutable CHNSLView chnsOldLocal;
+
+ Param& param; // parameter functions
+
+public:
+ enum {doPatternVolume=true};
+ //enum {doLambdaVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ THINSTRIP_FULL_CHNS_Monolithic(Param& param_, CHNSGFS& chnsGfs_, CHNSContainer& chnsVOld_) :
+ chnsOldLocal(chnsGfs_,chnsVOld_),
+ param(param_)
+ {
+ }
+
+
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+ {
+ //Quadrature Rule
+ const int dim = EG::Entity::dimension;
+ const int order = 2; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chnsOldLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHNSVariables_All<RF,dim,LFSV> chns(chCache, vCache, pCache, ip.position(), lfsv, x.base(), S);
+ CHNSVariables_All<RF,dim,LFSV> chnsOld(chCache, vCache, pCache, ip.position(), lfsv, chnsOldLocal.lv, S);
+
+ RF pf3 = 1-chns.pf-chns.pf2;
+ RF xi3 = -param.tw.Sigma3*(chns.xi/param.tw.Sigma1 + chns.xi2/param.tw.Sigma2);
+ RF pf_f= chns.pf+chns.pf2+2*param.delta*pf3;
+
+ RF pf_c = chns.pf;
+ RF pf_c_Old = chnsOld.pf;
+ RF pf_c_reg = pf_c+param.delta;
+ RF pf_c_reg_Old = pf_c_Old + param.delta;
+
+
+ Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ Dune::FieldVector<RF,dim> outerGradP(0.0);
+ outerGradP[0]=-param.gradP;
+ gradpf_f.axpy((1-2*param.delta),chns.gradpf+chns.gradpf2);
+
+ RF rho_f=param.phys.rho_f(chns.pf,chns.pf2);
+ RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
+
+ for (size_t i=0; i<chns.pfspace.size(); i++)
+ {
+ //RF Rspeed = std::max(param.tw.eval_q(ch.pf,ch.pf2,pf3) - 0.1, 0.0);
+ RF Rspeed = param.tw.eval_q(chns.pf,chns.pf2,pf3);
+ RF R1 = -Rspeed/param.eps*
+ (param.reaction.eval(chns.u)+param.phys.curvatureAlpha*(chns.xi-xi3))*chns.basis.phi[i];
+ RF J1 = param.mobility.eval(chns.pf,chns.pf2,pf3)*param.eps/param.tw.Sigma1 * (chns.gradxi*chns.basis.gradphi[i][0]);
+
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(chns.pfspace,i,factor*(
+ (chns.pf-chnsOld.pf)*chns.basis.phi[i]/param.time.dt
+ +J1
+ +(chns.pf*chns.div_v + (chns.gradpf*chns.v))*chns.basis.phi[i] //TODO div_v or div_v OLD?, v or vOld
+ //-pf1*(v*gradphi[i][0]) weak transport
+ -R1
+ ));
+
+ r.accumulate(chns.pf2space,i,factor*(
+ (chns.pf2-chnsOld.pf2)*chns.basis.phi[i]/param.time.dt
+ +param.mobility.eval(chns.pf,chns.pf2,pf3)*param.eps/param.tw.Sigma2 * (chns.gradxi2*chns.basis.gradphi[i][0])
+ +(chns.pf2*chns.div_v + (chns.gradpf2*chns.v))*chns.basis.phi[i] //TODO div_v or div_v OLD?, v or vOld
+ //-pf2*(v*gradphi[i][0])
+ ));
+
+ r.accumulate(chns.xispace,i,factor*(
+ chns.xi*chns.basis.phi[i]
+ - 1/param.eps* param.tw.eval_dpf1(chns.pf,chns.pf2,pf3) *chns.basis.phi[i]
+ - param.eps * param.tw.Sigma1* (chns.gradpf*chns.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(chns.xi2space,i,factor*(
+ chns.xi2*chns.basis.phi[i]
+ - 1/param.eps* param.tw.eval_dpf2(chns.pf,chns.pf2,pf3) *chns.basis.phi[i]
+ - param.eps * param.tw.Sigma2* (chns.gradpf2*chns.basis.gradphi[i][0])
+ ));
+
+ // integrate (pf-pfold)*(u-uast)*phi_i/delta_t + pfold*(u-uold)*phi_i/delta_t + D*pf*gradu*grad phi_i = 0
+
+// r.accumulate(chns.uspace,i,factor*(chns.u-param.phys.uAst)*chns.basis.phi[i]);
+
+
+
+ //Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
+ r.accumulate(chns.uspace,i,factor*(
+
+ //(u-param.uAst)*phi[i]/param.dt //For Testing
+
+ (pf_c_reg-pf_c_reg_Old)*(chns.u-param.phys.uAst)*chns.basis.phi[i]/param.time.dt
+ + pf_c_reg_Old*(chns.u-chnsOld.u)*chns.basis.phi[i]/param.time.dt
+ + param.phys.D*pf_c_reg_Old*(chns.gradu*chns.basis.gradphi[i][0])
+ + J1*(chns.u-param.phys.uAst)
+ //- pf_c*(u-param.phys.uAst)*(v*gradphi[i][0])
+ + (pf_c*(chns.u-param.phys.uAst)*chns.div_v //TODO v or vOld
+ + (chns.u-param.phys.uAst)*(chns.gradpf*chns.v)
+ + pf_c*(chns.gradu*chns.v) )*chns.basis.phi[i]
+ ));
+ } // for i
+
+/* for (size_t i=0; i<chns.pspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ r.accumulate(chns.pspace,i, -1.0 * (chns.p )* chns.pbasis.phi[i] * factor);
+ }
+ }
+ for (size_t i=0; i<chns.vfemspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ r.accumulate(child(chns.vspace,k),i, ( chns.v[k]*chns.vbasis.phi[i] )*factor);
+ } //for k
+ } // for i
+*/
+
+
+ for (size_t i=0; i<chns.pspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: (div v)*phi
+ //r.accumulate(ns.pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor); Weak form
+ r.accumulate(chns.pspace,i, -1.0 * (pf_f*chns.div_v + (gradpf_f*chns.v) )* chns.pbasis.phi[i] * factor);
+ }
+ }
+ for (size_t i=0; i<chns.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*(chns.gradxi[k]/param.tw.Sigma1 + chns.gradxi2[k]/param.tw.Sigma2);
+ RF S = 0;
+ S += - chns.xi2*(chns.gradpf[k]*pf_f - chns.pf * gradpf_f[k])/(pf_f*pf_f); //TODO!!!!!!! Other pf_f scaling!!!!!
+ S += - chns.xi *(chns.gradpf2[k]*pf_f - chns.pf2 * gradpf_f[k])/(pf_f*pf_f);
+ S += - 2*param.delta*pf3*(gradxi3k-chns.gradxi2[k]-chns.gradxi[k]);
+
+ //Stronger form - flow in interfaces possible
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k] - 2*param.delta*pf3*gradxi3k;
+
+ //Works but is not the one in the theory (same as above, w/o regularization)
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k];
+
+ S = S*param.SurfaceTensionScaling*chns.vbasis.phi[i];
+
+ r.accumulate(child(chns.vspace,k),i, ( rho_f_reg*(chns.v[k]-chnsOld.v[k])*chns.vbasis.phi[i]/param.time.dt //TimeEvolution
+ //- ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k]) //Pressure in weak form
+ + pf_f * ((chns.gradp[k]+outerGradP[k]) * chns.vbasis.phi[i]) //Pressure in strong form
+ + param.phys.eval_mu(pf_f) * (chns.jacv[k]*chns.vbasis.gradphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.vDis.eval(pf_f)*chns.v[k]*chns.vbasis.phi[i] //Velocity Dissipation in Solid
+ //- rho*g[k]*vphi[i] //Gravitation
+ - S //Surface Tension
+ )*factor);
+ if(k==0)
+ {
+ r.accumulate(child(chns.vspace,k),i, ( pf_f * (param.phys.outerMomentumXForce * chns.vbasis.phi[i]) // Momentum Forcing Term
+ )*factor);
+ }
+ } //for k
+ } // for i
+
+ } // for ip
+ }
+
+};