common_functionals.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/analytic/function.hpp>
#include <kernel/assembly/linear_functional.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    namespace Common
    {
      template<typename Function_>
      class ForceFunctional :
        public LinearFunctional
      {
      public:
        static_assert(Function_::can_value, "function can't compute values");
 
        static constexpr TrafoTags trafo_config = TrafoTags::img_point;
 
        static constexpr SpaceTags test_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public LinearFunctional::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::TrafoEvaluator TrafoEvaluator;
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename Analytic::EvalTraits<DataType, Function_>::ValueType ValueType;
 
        protected:
          typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
          typename Function_::template Evaluator<AnalyticEvalTraits> _func_eval;
          typename AnalyticEvalTraits::ValueType _func_value;
 
        public:
          explicit Evaluator(const ForceFunctional<Function_>& functional) :
            _func_eval(functional._function),
            _func_value(DataType(0))
          {
          }
 
          void set_point(const TrafoData& tau)
          {
            // evaluate function value
            _func_value = _func_eval.value(tau.img_point);
          }
 
          // copy pasted since Doxygen does not like the operator part in
          // \copydoc LinearFunctional::Evaluator::operator()
          ValueType eval(const TestBasisData& psi) const
          {
            return _func_value * psi.value;
          }
        }; // class ForceFunctional::Evaluator<...>
 
      protected:
        const Function_& _function;
 
      public:
        explicit ForceFunctional(const Function_& function) :
          _function(function)
        {
        }
      }; // class ForceFunctional
 
      template<typename Function_>
      class LaplaceFunctional :
        public LinearFunctional
      {
      public:
        static_assert(Function_::can_hess, "function can't compute hessians");
 
        static constexpr TrafoTags trafo_config = TrafoTags::img_point;
        static constexpr SpaceTags test_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public LinearFunctional::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::TrafoEvaluator TrafoEvaluator;
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename Analytic::EvalTraits<DataType, Function_>::ValueType ValueType;
 
        protected:
          typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
          typename Function_::template Evaluator<AnalyticEvalTraits> _func_eval;
          typename AnalyticEvalTraits::ValueType _func_laplace;
 
        public:
          explicit Evaluator(const LaplaceFunctional<Function_>& functional) :
            _func_eval(functional._function),
            _func_laplace(DataType(0))
          {
          }
 
          void set_point(const TrafoData& tau)
          {
            // compute the function hessian's trace (=laplacian)
            //_func_laplace = - _func_eval.hessian(tau.img_point).trace();
            calc_trace(_func_eval.hessian(tau.img_point));
          }
 
          template<typename T_, int m_, int n_, int sm_, int sn_>
          void calc_trace(const Tiny::Matrix<T_, m_, n_, sm_, sn_>& hess)
          {
            _func_laplace = - hess.trace();
          }
 
          template<typename T_, int l_, int m_, int n_, int sl_, int sm_, int sn_>
          void calc_trace(const Tiny::Tensor3<T_, l_, m_, n_, sl_, sm_, sn_>& hess)
          {
            _func_laplace.format();
            for(int i(0); i < hess.l; ++i)
              for(int j(0); j < hess.m; ++j)
                _func_laplace[i] -= hess[i][j][j];
          }
 
          // copy pasted since Doxygen does not like the operator part in
          // \copydoc LinearFunctional::Evaluator::operator()
          ValueType eval(const TestBasisData& psi) const
          {
            return _func_laplace * psi.value;
          }
        }; // class LaplaceFunctional::Evaluator<...>
 
      protected:
        const Function_& _function;
 
      public:
        explicit LaplaceFunctional(const Function_& function) :
          _function(function)
        {
        }
      }; // class LaplaceFunctional
    } // namespace Common
  } // namespace Assembly
} // namespace FEAT