Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
ldomain.cc 21.57 KiB
1 
// -*- tab-width: 2; indent-tabs-mode: nil -*-
2 
#ifdef HAVE_CONFIG_H
3 
#include "config.h"
4 
#endif
5
6 
#include <dune/pdelab/boilerplate/pdelab.hh>
7 
#include <dune/pdelab/localoperator/convectiondiffusionfem.hh>
8
9 
//***********************************************************************
10 
//***********************************************************************
11 
// define the reentrant corner in the L-domain with known exact solution
12 
//***********************************************************************
13 
//***********************************************************************
14
15 
template<typename GV, typename RF>
16 
class ReentrantCornerProblem
17 
{
18 
  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
19
20 
public:
21 
  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
22
23 
  //! tensor diffusion coefficient
24 
  typename Traits::PermTensorType
25 
  A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
26 
  {
27 
    typename Traits::PermTensorType I;
28 
    for (std::size_t i=0; i<Traits::dimDomain; i++)
29 
      for (std::size_t j=0; j<Traits::dimDomain; j++)
30 
        I[i][j] = (i==j) ? 1.0 : 0.0;
31 
    return I;
32 
  }
33
34 
  //! velocity field
35 
  typename Traits::RangeType
36 
  b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
37 
  {
38 
    typename Traits::RangeType v(0.0);
39 
    return v;
40 
  }
41
42 
  //! sink term
43 
  typename Traits::RangeFieldType
44 
  c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
45 
  {
46 
    return 0.0;
47 
  }
48
49 
  //! source term
50 
  typename Traits::RangeFieldType
51 
  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
52 
  {
53 
    return 0.0;
54 
  }
55
56 
  //! boundary condition type function
57 
  BCType
58 
  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
59 
  {
60 
    //typename Traits::DomainType xglobal = is.geometry().global(x);
61 
    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
62 
  }
63
64 
  //! Dirichlet boundary condition value
65 
  typename Traits::RangeFieldType
66 
  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
67 
  {
68 
    typename Traits::DomainType x = e.geometry().global(xlocal);
69
70 
    typename Traits::DomainFieldType theta = std::atan2(x[1], x[0]);
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
if(theta < 0.0) theta += 2*M_PI; typename Traits::DomainFieldType r = x.two_norm(); return pow(r,2.0/3.0)*std::sin(theta*2.0/3.0); } //! Neumann boundary condition typename Traits::RangeFieldType j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const { return 0.0; } //! Neumann boundary condition typename Traits::RangeFieldType o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const { return 0.0; } }; //*********************************************************************** //*********************************************************************** // a grid function giving the gradient of the exact solution // needed for computing H^1 errors //*********************************************************************** //*********************************************************************** //! exact gradient of solution template<typename GV, typename RF> class ExactGradient : public Dune::PDELab::GridFunctionBase< Dune::PDELab::GridFunctionTraits<GV,RF, GV::dimension,Dune::FieldVector<RF,GV::dimension> >, ExactGradient<GV,RF> > { GV gv; public: typedef Dune::PDELab::GridFunctionTraits<GV,RF, GV::dimension,Dune::FieldVector<RF,GV::dimension> > Traits; typedef Dune::PDELab::GridFunctionBase<Traits,ExactGradient<GV,RF> > BaseT; ExactGradient (const GV& gv_) : gv(gv_) {} inline void evaluate (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal, typename Traits::RangeType& y) const { typename Traits::DomainType x = e.geometry().global(xlocal); typename Traits::DomainFieldType theta = std::atan2(x[1], x[0]); if(theta < 0.0) theta += 2*M_PI; typename Traits::DomainFieldType r = x.two_norm(); y[0] = (2.0/3.0)*pow(r,-1.0/3.0)*(cos(theta)*sin(2.0*theta/3.0) - sin(theta)*cos(2.0*theta/3.0)); y[1] = (2.0/3.0)*pow(r,-1.0/3.0)*(sin(theta)*sin(2.0*theta/3.0) + cos(theta)*cos(2.0*theta/3.0)); } inline const typename Traits::GridViewType& getGridView () const { return gv; } }; //*********************************************************************** //*********************************************************************** // some function adapters //*********************************************************************** //***********************************************************************
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
/*! \brief Adapter returning f1(x)-f2(x) for two given grid functions \tparam T1 a grid function type \tparam T2 a grid function type */ template<typename T1, typename T2> class DifferenceAdapter : public Dune::PDELab::GridFunctionBase< Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType, typename T1::Traits::RangeFieldType, 1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > ,DifferenceAdapter<T1,T2> > { public: typedef Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType, typename T1::Traits::RangeFieldType, 1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > Traits; //! constructor DifferenceAdapter (const T1& t1_, const T2& t2_) : t1(t1_), t2(t2_) {} //! \copydoc GridFunctionBase::evaluate() inline void evaluate (const typename Traits::ElementType& e, const typename Traits::DomainType& x, typename Traits::RangeType& y) const { typename Traits::RangeType y1; t1.evaluate(e,x,y1); typename Traits::RangeType y2; t2.evaluate(e,x,y2); y1 -= y2; y = y1; } inline const typename Traits::GridViewType& getGridView () const { return t1.getGridView(); } private: const T1& t1; const T2& t2; }; /*! \brief Adapter returning ||f1(x)-f2(x)||^2 for two given grid functions \tparam T1 a grid function type \tparam T2 a grid function type */ template<typename T1, typename T2> class DifferenceSquaredAdapter : public Dune::PDELab::GridFunctionBase< Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType, typename T1::Traits::RangeFieldType, 1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > ,DifferenceSquaredAdapter<T1,T2> > { public: typedef Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType, typename T1::Traits::RangeFieldType, 1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > Traits; //! constructor DifferenceSquaredAdapter (const T1& t1_, const T2& t2_) : t1(t1_), t2(t2_) {} //! \copydoc GridFunctionBase::evaluate() inline void evaluate (const typename Traits::ElementType& e, const typename Traits::DomainType& x, typename Traits::RangeType& y) const