laplace.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) 2012-2019 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef LAPLACE_OP_HH
#define LAPLACE_OP_HH
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
#include "fem/fixdune.hh"
#include "utilities/linalg/scalarproducts.hh"
 
namespace Kaskade 
{
 
  template <class Scalar, int dim, bool isotropic=true>
  class ScalarLaplace 
  {
  public:
    using DiffusionTensor = std::conditional_t<isotropic,
                                               Scalar,
                                               Dune::FieldMatrix<Scalar,dim,dim>>;
    ScalarLaplace()
    {
      if constexpr (isotropic)
        kappa = 1;
      else
        kappa = unitMatrix<Scalar,dim,dim>();
    }
    
    Laplace& setDiffusionTensor(DiffusionTensor const& kappa_)
    {
      kappa = kappa_;
      return *this;
    }
    
    Laplace& setLinearizationPoint(Dune::FieldVector<Scalar,dim> const& du0_)
    {
      du0 = du0_;
      return *this;
    }
    
    Scalar d0() const
    {
      auto a = kappa*du0;
      return 0.5*(a*du0);
    }
    
    Dune::FieldVector<Scalar,1> d1(VariationalArg<Scalar,dim> const& arg) const
    {
      auto a = kappa*du0;
      return arg.derivative*a;
    }
    
    Dune::FieldMatrix<Scalar,1,1> d2(VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const 
    {
      auto a = kappa*arg1.derivative[0];
      return a*arg2.derivative[0];
    }
    
  private:
    DiffusionTensor kappa;               // diffusion constant
    Dune::FieldVector<Scalar,dim> du0;   // gradient of linearization point
  };
 
  template <class Scalar, int dim, int components, bool isotropic=true>
  class Laplace 
  {
  public:
    Laplace(): kappa(1), du0(0)
    {
    }
 
    void setDiffusionTensor(Scalar kappa_)
    {
      // TODO: shall we enforce kappa > 0? And kappa real?
      kappa = kappa_;
    }
 
    void setLinearizationPoint(Dune::FieldMatrix<Scalar,components,dim> const& du0_)
    {
      du0 = du0_;
    }
 
    Scalar d0() const
    {
      return 0.5*kappa*sp(du0, du0);
    }
 
    Scalar d1(VariationalArg<Scalar,dim,components> const& arg) const
    {
      return kappa*sp(arg.derivative, du0);
    }
 
    Scalar d2(VariationalArg<Scalar,dim,components> const &arg1, VariationalArg<Scalar,dim,components> const &arg2) const
    {
      return kappa*sp(arg1.derivative, arg2.derivative);
    }
 
  private:
    // diffusion constant
    Scalar kappa;
    Dune::FieldMatrix<Scalar,components,dim> du0;
    LinAlg::EuclideanScalarProduct sp;
  };
 
 
 
}
 
#endif