diff --git a/dune/phasefield/fem_2p_cahnhilliard_a.hh b/dune/phasefield/fem_2p_cahnhilliard_a.hh new file mode 100644 index 0000000000000000000000000000000000000000..b851991755db4c3fae915659c8fcaee538d60330 --- /dev/null +++ b/dune/phasefield/fem_2p_cahnhilliard_a.hh @@ -0,0 +1,89 @@ +#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_2P_CahnHilliard_a : + public Dune::PDELab::FullVolumePattern, + public Dune::PDELab::LocalOperatorDefaultFlags, + public Dune::PDELab::NumericalJacobianVolume<FEM_2P_CahnHilliard_a<Param, CHGFS, CHContainer> >, + public Dune::PDELab::NumericalJacobianApplyVolume<FEM_2P_CahnHilliard_a<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_2P_CahnHilliard_a(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) : + chOldLocal(chGfs_,chVOld_), + param(param_) + { + } + + + template<typename EG, typename LFSU, typename X, typename LFSV, typename R> + void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const + { + //Quadrature Rule + const int dim = EG::Entity::dimension; + const int order = 5; + 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 J1 = param.mobility.eval(ch.pf) * (ch.gradxi*ch.basis.gradphi[i][0]); + + const float alpha = 50; + + // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ... + //weak form CH pf=c, xi=z, phi=testfkt + r.accumulate(ch.pfspace,i,factor*( + (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt + + 1/param.eps* param.dw.eval_dd(ch.pf) * (ch.gadpf*ch.basis.gradphi[i][0]) + - alpha * (ch.gradxi*ch.basis.gradphi[i][0] ) + + alpha * (ch.gradpf*ch.basis.gradphi[i][0]) + )); +//weak from laplace z + r.accumulate(ch.xispace,i,factor*( + ch.gradxi*ch.basis.gradphi[i][0] + + alpha (ch.xi*ch.basis.phi[i] ) + - alpha (ch.pf*ch.basis.phi[i] ) + )); + } + } + } + +};