mmdgproblem5.hh

#ifndef __DUNE_MMDG_MMDGPROBLEM5_HH__
#define __DUNE_MMDG_MMDGPROBLEM5_HH__

#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>

//a coupled dg problem
template<class Coordinate, class Scalar = double>
class MMDGProblem5 : public CoupledDGProblem<Coordinate, Scalar>
{
public:
  static constexpr int dim = Coordinate::dimension;
  using Base = CoupledDGProblem<Coordinate, Scalar>;
  using Matrix = typename Base::Matrix;

  //constructor
  MMDGProblem5 (const Scalar d) : d_(d) {}

  //the exact bulk solution at position pos
  Scalar exactSolution (const Coordinate& pos) const
  {
    return (pos[1] > pos[0]) ?
      1.0 / (1.0 + pos[0] * pos[0]) : -1.0 / (1.0 + pos[1] * pos[1]);
  }

  //the exact solution on the interface at position pos
  Scalar exactInterfaceSolution (const Coordinate& pos) const
  {
    return 0.0;
  }

  //indicates whether an exact solution is implemented for the problem
  bool hasExactSolution () const
  {
    return true;
  };

  //source term at position pos
  Scalar q (const Coordinate& pos) const
  {
    const Scalar x2Plus1 = pos[0] * pos[0] + 1.0;
    const Scalar y2Plus1 = pos[1] * pos[1] + 1.0;

    return (pos[1] > pos[0]) ?
      -4.0 * pos[0] / (x2Plus1 * x2Plus1) : 4.0 * pos[1] / (y2Plus1 * y2Plus1);
  }

  //interface source term at position pos
  Scalar qInterface (const Coordinate& pos) const
  {
    return 0.0;
  }

  //aperture d of the fracture at position pos
  Scalar aperture (const Coordinate& pos) const
  {
    return d_;
  }

  //permeability tensor at position pos
  Matrix K (const Coordinate& pos) const
  {
    Matrix permeability(0.0);

    for (int i = 0; i < dim; i++)
    {
      permeability[i][i] = pos[i] + 1.0 / pos[i];
    }

    return permeability;
  }

  //permeability of the fracture in normal direction at position pos
  Scalar Kperp (const Coordinate& pos) const
  {
    constexpr Scalar Kvalue = 1.0 / sqrt(2.0);
    return Kvalue;
  }

  //tangential permeability tensor of the interface at position pos
  Matrix Kparallel (const Coordinate& pos) const
  {
    Matrix permeability(0.0);
    Scalar y2Plus1_2 = 1.0 + pos[1] * pos[1];
    y2Plus1_2 *= y2Plus1_2;

    for (int i = 0; i < dim; i++)
    {
      permeability[i][i] = y2Plus1_2;
    }

    return permeability;
  }

  //returns the recommended quadrature order to compute an integral
  //over x * boundary(x)
  int quadratureOrderBoundary () const
  {
    return 10;
  }

  //returns the recommended quadrature order to compute an integral
  //over x * interfaceBoundary(x)
  int quadratureOrderInterfaceBoundary () const
  { //NOTE: only relevant for 3d, not implemented!
    return 10;
  }

  //returns the recommended quadrature order to compute an integral
  //over x * qInterface(x)
  int quadratureOrder_qInterface () const
  {
    return 10;
  }

  //returns the recommended quadrature order to compute an integral
  //over x * q(x)
  int quadratureOrder_q () const
  {
    return 10;
  }

private:
  Scalar d_; //aperture

};

#endif