evaluator_base.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/trafo/eval_data.hpp>
 
namespace FEAT
{
  namespace Trafo
  {
    namespace Intern
    {
      template<bool _enable>
      struct TrafoEvalHelper;
    } // namespace Intern
 
    template<
      typename Trafo_,
      typename Evaluator_,
      typename EvalPolicy_>
    class EvaluatorBase
    {
    public:
      typedef Trafo_ TrafoType;
      typedef typename TrafoType::ShapeType ShapeType;
 
      typedef EvalPolicy_ EvalPolicy;
 
      typedef EvalPolicy_ EvalTraits;
 
      typedef typename EvalPolicy::DataType DataType;
      typedef typename EvalPolicy::DomainPointType DomainPointType;
      typedef typename EvalPolicy::ImagePointType ImagePointType;
      typedef typename EvalPolicy::JacobianMatrixType JacobianMatrixType;
      typedef typename EvalPolicy::JacobianInverseType JacobianInverseType;
      typedef typename EvalPolicy::JacobianDeterminantType JacobianDeterminantType;
      typedef typename EvalPolicy::HessianTensorType HessianTensorType;
      typedef typename EvalPolicy::HessianInverseType HessianInverseType;
 
      static constexpr int domain_dim = EvalPolicy::domain_dim;
      static constexpr int image_dim = EvalPolicy::image_dim;
 
      typedef Index CellIterator;
 
      // Note:
      // The following block serves as an element interface documentation and is therefore only
      // visible to doxygen. The actual functionality has to be supplied by the implementation.
#ifdef DOXYGEN
      static constexpr TrafoTags eval_caps = ...;
#endif // DOXYGEN
 
      template<TrafoTags cfg_>
      struct ConfigTraits
      {
        static constexpr TrafoTags config =
          (*(cfg_ &  TrafoTags::hess_inv) ?
            TrafoTags::hess_inv : TrafoTags::none) |
          (*(cfg_ & (TrafoTags::hess_ten | TrafoTags::hess_inv)) ?
            TrafoTags::hess_ten : TrafoTags::none) |
          (*(cfg_ & TrafoTags::jac_det) ?
            TrafoTags::jac_det : TrafoTags::none) |
          (*(cfg_ & (TrafoTags::jac_inv | TrafoTags::hess_inv)) ?
            TrafoTags::jac_inv : TrafoTags::none) |
          //(*(cfg_ & (TrafoTags::jac_mat | TrafoTags::jac_det | TrafoTags::jac_inv | TrafoTags::hess_inv | TrafoTags::normal)) ?
          (*(cfg_ & (TrafoTags::jac_mat | TrafoTags::jac_det | TrafoTags::jac_inv | TrafoTags::hess_inv)) ?
            TrafoTags::jac_mat : TrafoTags::none) |
          (*(cfg_ & TrafoTags::normal) ?
            TrafoTags::normal : TrafoTags::none) |
          (*(cfg_ & (TrafoTags::img_point)) ?
            TrafoTags::img_point : TrafoTags::none) |
          TrafoTags::dom_point;
 
        typedef Trafo::EvalData<EvalTraits, config> EvalDataType;
      }; // struct ConfigTraits<...>
 
    protected:
      Evaluator_& cast()
      {
        return static_cast<Evaluator_&>(*this);
      }
 
      const Evaluator_& cast() const
      {
        return static_cast<const Evaluator_&>(*this);
      }
 
      const TrafoType& _trafo;
 
      Index _cell_index;
 
      explicit EvaluatorBase(const TrafoType& trafo) :
        _trafo(trafo),
        _cell_index(~Index(0))
      {
      }
 
    public:
      const TrafoType& get_trafo() const
      {
        return _trafo;
      }
 
      CellIterator begin() const
      {
        return 0;
      }
 
      CellIterator end() const
      {
        return get_num_cells();
      }
 
      Index get_num_cells() const
      {
        return _trafo.get_mesh().get_num_entities(domain_dim);
      }
 
      Index get_cell_index() const
      {
        return _cell_index;
      }
 
      void prepare(const CellIterator& cell)
      {
        // store cell index
        _cell_index = cell;
      }
 
      void finish()
      {
        // reset cell index
        _cell_index = ~Index(0);
      }
 
      template<TrafoTags cfg_>
      void operator()(Trafo::EvalData<EvalTraits, cfg_>& trafo_data, const DomainPointType& dom_point) const
      {
        // typedef mumbo-jumbo
        typedef Trafo::EvalData<EvalTraits, cfg_> TrafoData;
 
        // Note:
        // The following static constexpr bool variables are required to
        // circumvent a compiler bug in GCC 6.1
 
        // store domain point
        static constexpr bool want_dom_point = *(TrafoData::config & TrafoTags::dom_point);
        Intern::TrafoEvalHelper<want_dom_point>::set_dom_point(trafo_data, dom_point);
        // map image point
        static constexpr bool want_img_point = *(TrafoData::config & TrafoTags::img_point);
        Intern::TrafoEvalHelper<want_img_point>::map_img_point(trafo_data, cast());
        // calculate jacobian matrix
        static constexpr bool want_jac_mat = *(TrafoData::config & TrafoTags::jac_mat);
        Intern::TrafoEvalHelper<want_jac_mat>::calc_jac_mat(trafo_data, cast());
        // calculate inverse jacobian matrix
        static constexpr bool want_jac_inv = *(TrafoData::config & TrafoTags::jac_inv);
        Intern::TrafoEvalHelper<want_jac_inv>::calc_jac_inv(trafo_data, cast());
        // calculate jacobian determinants
        static constexpr bool want_jac_det = *(TrafoData::config & TrafoTags::jac_det);
        Intern::TrafoEvalHelper<want_jac_det>::calc_jac_det(trafo_data, cast());
        // calculate hessian tensor
        static constexpr bool want_hess_ten = *(TrafoData::config & TrafoTags::hess_ten);
        Intern::TrafoEvalHelper<want_hess_ten>::calc_hess_ten(trafo_data, cast());
        // calculate inverse hessian tensor
        static constexpr bool want_hess_inv = *(TrafoData::config & TrafoTags::hess_inv);
        Intern::TrafoEvalHelper<want_hess_inv>::calc_hess_inv(trafo_data, cast());
 
        // note: normal vectors are computed directly by the assembler classes
        // calculate normal vector
        //static constexpr bool want_normal = *(TrafoData::config & TrafoTags::normal);
        //Intern::TrafoEvalHelper<want_normal>::calc_normal(trafo_data, cast());
      }
 
      // Note:
      // The following block serves as an element interface documentation and is therefore only
      // visible to doxygen. The actual functionality has to be supplied by the implementation.
#ifdef DOXYGEN
      void map_point(ImagePointType& img_point, const DomainPointType& dom_point) const;
 
      void calc_jac_mat(JacobianMatrixType& jac_mat, const DomainPointType& dom_point) const;
 
      void calc_hess_ten(HessianTensorType& hess_ten, const DomainPointType& dom_point) const;
 
      DataType volume() const;
 
      DataType width_directed(const ImagePointType& ray) const;
#endif // DOXYGEN
    }; // class EvaluatorBase<...>
 
    namespace Intern
    {
      template<bool _enable>
      struct TrafoEvalHelper
      {
        template<typename TrafoData_, typename DomPoint_>
        static void set_dom_point(TrafoData_&, const DomPoint_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void map_img_point(TrafoData_&, const Evaluator_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_mat(TrafoData_&, const Evaluator_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_inv(TrafoData_&, const Evaluator_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_det(TrafoData_&, const Evaluator_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_hess_ten(TrafoData_&, const Evaluator_&) {}
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_hess_inv(TrafoData_&, const Evaluator_&) {}
 
        //template<typename TrafoData_, typename Evaluator_>
        //static void calc_normal(TrafoData_&, const Evaluator_&) {}
      };
 
      template<>
      struct TrafoEvalHelper<true>
      {
        template<typename TrafoData_, typename DomPoint_>
        static void set_dom_point(TrafoData_& trafo_data, const DomPoint_& dom_point)
        {
          trafo_data.dom_point = dom_point;
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void map_img_point(TrafoData_& trafo_data, const Evaluator_& evaluator)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::img_point))
            XABORTM("trafo evaluator can't compute image point coordinates");
          evaluator.map_point(trafo_data.img_point, trafo_data.dom_point);
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_mat(TrafoData_& trafo_data, const Evaluator_& evaluator)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::jac_mat))
            XABORTM("trafo evaluator can't compute jacobian matrices");
          // let the evaluator compute the jacobian matrix
          evaluator.calc_jac_mat(trafo_data.jac_mat, trafo_data.dom_point);
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_inv(TrafoData_& trafo_data, const Evaluator_&)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::jac_inv))
            XABORTM("trafo evaluator can't compute jacobian inverse matrices");
          // invert the jacobian matrix
          trafo_data.jac_inv.set_inverse(trafo_data.jac_mat);
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_jac_det(TrafoData_& trafo_data, const Evaluator_&)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::jac_det))
            XABORTM("trafo evaluator can't compute jacobian determinants");
          // compute the volume of the jacobian matrix
          trafo_data.jac_det = trafo_data.jac_mat.vol();
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_hess_ten(TrafoData_& trafo_data, const Evaluator_& evaluator)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::hess_ten))
            XABORTM("trafo evaluator can't compute hessian tensors");
          // let the evaluator compute the hessian tensor
          evaluator.calc_hess_ten(trafo_data.hess_ten, trafo_data.dom_point);
        }
 
        template<typename TrafoData_, typename Evaluator_>
        static void calc_hess_inv(TrafoData_& trafo_data, const Evaluator_&)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::hess_inv))
            XABORTM("trafo evaluator can't compute inverse hessian tensors");
 
          typedef typename TrafoData_::EvalTraits EvalTraits;
          typedef typename EvalTraits::DataType DataType;
          typename EvalTraits::HessianTensorType hess_jac;
 
          // compute inverse
          hess_jac.format();
          hess_jac.add_double_mat_mult(trafo_data.hess_ten, trafo_data.jac_inv, trafo_data.jac_inv);
          trafo_data.hess_inv.format();
          trafo_data.hess_inv.add_mat_tensor_mult(trafo_data.jac_inv, hess_jac, -DataType(1));
        }
 
        /*template<typename TrafoData_, typename Evaluator_>
        static void calc_normal(TrafoData_& trafo_data, const Evaluator_&)
        {
          if(!*(Evaluator_::eval_caps & TrafoTags::normal))
            XABORTM("trafo evaluator can't compute normal vectors");
          // compute the normal vector from of the jacobian matrix
          trafo_data.normal = Tiny::orthogonal(trafo_data.jac_mat);
          trafo_data.normal.normalize();
        }*/
      };
    } // namespace Intern
  } // namespace Trafo
} // namespace FEAT