boundaryConditions.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2019-2022 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef BOUNDARYCONDITIONS_HH
#define BOUNDARYCONDITIONS_HH
 
#include "fem/shapefunctioncache.hh"
 
namespace Kaskade
{
  template <class Functional>
  class HomogeneousNeumannBoundary
  {
    using Scalar = typename Functional::Scalar;
    static constexpr int dim = Functional::AnsatzVars::Grid::dimension;
 
  public:
    HomogeneousNeumannBoundary(Functional const&, typename Functional::OriginVars::VariableSet const&, int) {}
 
    HomogeneousNeumannBoundary(int) {}
 
    template <class FaceIterator>
    void moveTo(FaceIterator const&) {}
 
    template <class Evaluators>
    void evaluateAt(Dune::FieldVector<typename Functional::AnsatzVars::Grid::ctype, dim-1> const&, Evaluators const&) {}
 
    Scalar d0() const { return 0; }
 
    template <int row>
    auto d1(VariationalArg<Scalar,dim> const& ) const
    {
      return Dune::FieldVector<Scalar,Functional::TestVars::template Components<row>::m>(0);
    }
 
    template <int row, int col>
    auto d2(VariationalArg<Scalar,dim> const& , VariationalArg<Scalar,dim> const& ) const
    {
      return Dune::FieldMatrix<Scalar,Functional::TestVars::template Components<row>::m,Functional::AnsatzVars::template Components<row>::m>(0);
    }
  };
 
  template <class GridView, int components, class ScalarType=typename GridView::ctype>
  class DirichletPenaltyBoundary
  {
    using Scalar = ScalarType;
    static constexpr int dim = GridView::dimension;
    static constexpr int dimworld = GridView::dimensionworld;
 
  public:
 
    using Vector = Dune::FieldVector<Scalar,components>;
 
    void moveTo(typename GridView::IntersectionIterator const& fi)
    {
      // Compute h^{-2}. h is volume^(1/(dim-1)). Take care of special case for dim=1.
      if (dim > 1)
        weight = std::pow(fi->geometry().volume(),-2/(dim-1));
      else
        weight = 1;
    }
 
    void setBoundaryData(Scalar gamma_, Vector const& u_, Vector const& ud)
    {
      gamma = gamma_;
      assert(gamma>=0);
      u = u_;
      deviation = u-ud;
    }
 
    Scalar d0() const { return 0.5*gamma*weight * (deviation*deviation); }
 
    Vector d1(VariationalArg<Scalar,dim> const& argT) const
    {
      return gamma*weight*argT.value[0]*deviation;
    }
 
    Dune::FieldMatrix<Scalar,components,components>
    d2(VariationalArg<Scalar,dim> const& argTest, VariationalArg<Scalar,dim> const& argAnsatz) const
    {
      return gamma*weight*argTest.value[0]*argAnsatz.value[0] * unitMatrix<Scalar,components>();
    }
 
  private:
    Scalar gamma;
    Scalar weight; // this is h^{-2}
    Vector u, deviation;
  };
 
  template <class GridView, int components, class ScalarType=typename GridView::ctype>
  class DirichletNitscheBoundary
  {
    using Scalar = ScalarType;
    static constexpr int dim = GridView::dimension;
    static constexpr int dimworld = GridView::dimensionworld;
 
  public:
 
    using Vector = Dune::FieldVector<Scalar,components>;
 
    void moveTo(typename GridView::IntersectionIterator const& fi)
    {
      face = *fi;
 
      // Compute h^{-1}. h is volume^(1/(dim-1)). Take care of special case for dim=1.
      if (dim > 1)
        weight = std::pow(face.geometry().volume(),-1/(dim-1));
      else
        weight = 1;
    }
 
    void setBoundaryData(Dune::FieldVector<typename GridView::ctype, dim-1> const& xi,
                         Scalar gamma_, Vector const& u, Vector const& ud, Dune::FieldMatrix<Scalar,components,dim> ux_,
                         Dune::FieldMatrix<Scalar,dim,dim> const& kappa)
    {
      gamma = gamma_;
      assert(gamma>=0);
      deviation = u-ud;
      ux = ux_;
      kn = kappa * face.unitOuterNormal(xi);
    }
 
    Scalar d0() const { return 0.5*weight*gamma*(deviation*deviation) - deviation*(ux*kn); }
 
    Vector d1(VariationalArg<Scalar,dim> const& arg) const
    {
      return weight*gamma*arg.value[0]*deviation - arg.value[0]*(ux*kn) - deviation*(arg.derivative[0]*kn);
    }
 
    Dune::FieldMatrix<Scalar,components,components>
    d2(VariationalArg<Scalar,dim> const& argTest, VariationalArg<Scalar,dim> const& argAnsatz) const
    {
      Dune::FieldMatrix<Scalar,components,components> ones(1.0);
 
      Scalar t1 = argTest.value[0]*(argAnsatz.derivative[0]*kn);
      Scalar t2 = argAnsatz.value[0]*(argTest.derivative[0]*kn);
 
      return weight*gamma*argTest.value[0]*argAnsatz.value[0] * unitMatrix<Scalar,components>()
             - (t1+t2)*ones;
    }
 
  private:
    Scalar gamma;
    Scalar weight; // this is h^{-1}
    typename GridView::Intersection face;
    Dune::FieldVector<Scalar,dim> kn; // this is kappa*n
    Vector deviation;
    Dune::FieldMatrix<Scalar,components,dim> ux;
  };
 
 
}
 
#endif