Upload New File

dune/phasefield/fem_2p_cahnhilliard_a.hh

+#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] )
+        ));
+      }
+    }
+  }
+
+};