inverse_mapping.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
 
#include <kernel/trafo/mapping_base.hpp>
#include <kernel/util/exception.hpp>
 
#include <vector>
 
namespace FEAT
{
  namespace Trafo
  {
    namespace Intern
    {
      template<typename Shape_>
      struct InverseMappingHelper;
    } // namespace Intern
 
    class InverseMappingError :
      public FEAT::Exception
    {
    public:
      template<typename DT_, int n_, int s_>
      explicit InverseMappingError(const Tiny::Vector<DT_, n_, s_>& point, Index cell) :
        Exception(String("Failed to unmap point ") + stringify(point) + " on cell " + stringify(cell))
      {
      }
    }; // class InverseMappingError
 
    template<typename DataType_, int shape_dim_, int world_dim_ = shape_dim_>
    class InverseMappingData
    {
    public:
      typedef Tiny::Vector<DataType_, shape_dim_> DomainPointType;
      typedef Tiny::Vector<DataType_, world_dim_> ImagePointType;
 
      ImagePointType img_point;
 
      std::vector<Index> cells;
 
      std::vector<DomainPointType> dom_points;
 
      InverseMappingData()
      {
      }
 
      InverseMappingData(const InverseMappingData& other) = default;
 
      InverseMappingData(InverseMappingData&& other) :
        img_point(other.img_point),
        cells(std::forward<std::vector<Index>>(other.cells)),
        dom_points(std::forward<std::vector<DomainPointType>>(other.dom_points))
      {
      }
 
      InverseMappingData& operator=(InverseMappingData&& other)
      {
        if(this != &other)
        {
          img_point = other.img_point;
          cells = std::forward<std::vector<Index>>(other.cells);
          dom_points = std::forward<std::vector<DomainPointType>>(other.dom_points);
        }
        return *this;
      }
 
      bool empty() const
      {
        return cells.empty();
      }
 
      std::size_t size() const
      {
        return cells.size();
      }
    }; // class InverseMappingData
 
    template<typename Trafo_, typename DataType_>
    class InverseMapping
    {
    public:
      typedef Trafo_ TrafoType;
      typedef DataType_ DataType;
      typedef typename TrafoType::ShapeType ShapeType;
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int world_dim = TrafoType::world_dim;
 
      typedef InverseMappingData<DataType, shape_dim, world_dim> InvMapDataType;
 
      typedef typename InvMapDataType::ImagePointType ImagePointType;
      typedef typename InvMapDataType::DomainPointType DomainPointType;
 
    protected:
      typedef Tiny::Matrix<DataType_, 2, world_dim> BBox;
      const TrafoType& _trafo;
      std::vector<BBox> _bboxes;
      DataType _domain_tol;
      DataType _newton_tol;
      Index _newton_max_iter;
 
    public:
      explicit InverseMapping(const TrafoType& trafo,
        DataType bbox_tol = DataType(1E-2),
        DataType domain_tol = DataType(1E-4),
        bool init_boxes = true
        ) :
        _trafo(trafo),
        _domain_tol(domain_tol),
        _newton_tol(Math::pow(Math::eps<DataType>(), DataType(0.9))),
        _newton_max_iter(10u)
      {
        if(init_boxes)
          init_bounding_boxes(bbox_tol);
      }
 
      virtual ~InverseMapping()
      {
      }
 
      void init_bounding_boxes(DataType bbox_tol)
      {
        // get the mesh
        const auto& mesh = _trafo.get_mesh();
        // get the vertex set
        const auto& vtx_set = mesh.get_vertex_set();
        // get the vertices-at-cell index set
        const auto& vert_idx = mesh.template get_index_set<shape_dim, 0>();
 
        // resize bounding boxes array
        _bboxes.resize(std::size_t(mesh.get_num_elements()));
 
        // loop over all cells
        for(Index cell(0); cell < vert_idx.get_num_entities(); ++cell)
        {
          // get this cell's vertex indices
          const auto& vidx = vert_idx[cell];
 
          // get this cell's bounding box
          BBox& bbox = _bboxes.at(std::size_t(cell));
 
          // build cell bounding box
          bbox[0] = bbox[1] = vtx_set[vidx[0]];
          for(int i(1); i < vert_idx.get_num_indices(); ++i)
          {
            // get vertex and update bounding box
            const auto& vtx = vtx_set[vidx[i]];
            for(int j(0); j < world_dim; ++j)
            {
              Math::minimax(DataType(vtx[j]), bbox[0][j], bbox[1][j]);
            }
          }
 
          // finally, incorporate tolerances
          for(int j(0); j < world_dim; ++j)
          {
            // compute bbox dimension
            DataType bb_extra = bbox_tol * (bbox[1][j] - bbox[0][j]);
            bbox[0][j] -= bb_extra;
            bbox[1][j] += bb_extra;
          }
        }
      }
 
      InvMapDataType unmap_point(const ImagePointType& img_point, bool ignore_failures = false) const
      {
        // store image point
        InvMapDataType inv_data;
        inv_data.img_point = img_point;
 
        // find candidate cells
        std::vector<Index> cells;
        if(!this->find_candidate_cells(cells, img_point))
          return inv_data;
 
        // loop over all candidate cells
        for(Index cell : cells)
        {
          // try to unmap the point by newton
          DomainPointType dom_point;
          if(!this->unmap_point_by_newton(dom_point, img_point, cell))
          {
            // point unmapping failed!
            if(ignore_failures)
              continue;
 
            // throw an inverse mapping error
            throw InverseMappingError(dom_point, cell);
          }
 
          // check if the domain point is on the reference element
          if(!this->test_domain_point(dom_point))
          {
            // nope ==> skip point
            continue;
          }
 
          // okay
          inv_data.cells.push_back(cell);
          inv_data.dom_points.push_back(dom_point);
        }
 
        // found at least one cell?
        return inv_data;
      }
 
      InvMapDataType unmap_point(const ImagePointType& img_point, const std::vector<Index>& cells, bool ignore_failures = false) const
      {
        // store image point
        InvMapDataType inv_data;
        inv_data.img_point = img_point;
 
        if(cells.empty())
          return inv_data;
 
        // loop over all candidate cells
        for(Index cell : cells)
        {
          // try to unmap the point by newton
          DomainPointType dom_point;
          if(!this->unmap_point_by_newton(dom_point, img_point, cell))
          {
            // point unmapping failed!
            if(ignore_failures)
              continue;
 
            // throw an inverse mapping error
            throw InverseMappingError(dom_point, cell);
          }
 
          // check if the domain point is on the reference element
          if(!this->test_domain_point(dom_point))
          {
            // nope ==> skip point
            continue;
          }
 
          // okay
          inv_data.cells.push_back(cell);
          inv_data.dom_points.push_back(dom_point);
        }
 
        // found at least one cell?
        return inv_data;
      }
 
      bool find_candidate_cells(std::vector<Index>& cells, const ImagePointType& img_point) const
      {
        XASSERTM(_bboxes.size() == _trafo.get_mesh().get_num_elements(), "Boundingbox are not initilialized correctly");
        // loop over all cells/bounding boxes
        for(std::size_t cell(0); cell < _bboxes.size(); ++cell)
        {
          const BBox& bbox = _bboxes.at(cell);
 
          // check whether the point is in the bounding box
          bool is_in = true;
          for(int j(0); j < shape_dim; ++j)
          {
            if((img_point[j] < bbox[0][j]) || (img_point[j] > bbox[1][j]))
            {
              is_in = false;
              break;
            }
          }
 
          // point inside bounding box?
          if(is_in)
            cells.push_back(Index(cell));
        }
 
        return !cells.empty();
      }
 
      bool unmap_point_by_newton(DomainPointType& dom_point, const ImagePointType& img_point, const Index cell) const
      {
        // define an evaluator
        typedef typename TrafoType::template Evaluator<ShapeType, DataType>::Type TrafoEvaluator;
        TrafoEvaluator trafo_eval(_trafo);
 
        // for Newton iteration we need image points and inverse jacobians
        static constexpr TrafoTags trafo_tags = TrafoTags::img_point|TrafoTags::jac_inv;
 
        // define the trafo eval data
        typename TrafoEvaluator::template ConfigTraits<trafo_tags>::EvalDataType trafo_data;
 
        // initialize output domain point to cell centre
        for(int i(0); i < shape_dim; ++i)
          dom_point[i] = Shape::ReferenceCell<ShapeType>::template centre<DataType>(i);
 
        // compute squared tolerance
        const DataType tol_sqr = Math::sqr(_newton_tol);
 
        // no convergence yet
        bool converged = false;
 
        // prepare trafo evaluator
        trafo_eval.prepare(cell);
 
        // newton iteration
        for(Index iter(0); iter < _newton_max_iter; ++iter)
        {
          // evaluate trafo in current point
          trafo_eval(trafo_data, dom_point);
 
          // compute image point distance vector
          ImagePointType def = trafo_data.img_point - img_point;
 
          // check defect for tolerance
          if(def.norm_euclid_sqr() < tol_sqr)
          {
            converged = true;
            break;
          }
 
          // update reference point by newton
          dom_point -= trafo_data.jac_inv * def;
        }
 
        // finish trafo evaluator
        trafo_eval.finish();
 
        // Newton converged?
        return converged;
      }
 
      bool test_domain_point(const DomainPointType& dom_point) const
      {
        return Intern::InverseMappingHelper<ShapeType>::is_on_ref(dom_point, _domain_tol);
      }
    }; // class InverseMapping
 
    namespace Intern
    {
      template<int dim_>
      struct InverseMappingHelper<Shape::Simplex<dim_>>
      {
        // checks whether p is on the reference cell
        template<typename PT_, typename DT_>
        static bool is_on_ref(const PT_& p, const DT_ tol)
        {
          // reference simplex: all coords are in range [0,1]
          //                    and the sum of all coords is <= 1
          DT_ s = DT_(0);
          for(int i(0); i < dim_; ++i)
          {
            // coord outside range [-tol, 1+tol] ?
            if((p[i] < -tol) || (p[i] > DT_(1)+tol))
              return false;
            s += p[i];
          }
          // coord sum <= 1+tol ?
          return (s <= DT_(1) + tol);
        }
      };
 
      template<int dim_>
      struct InverseMappingHelper<Shape::Hypercube<dim_>>
      {
        // checks whether p is on the reference cell
        template<typename PT_, typename DT_>
        static bool is_on_ref(const PT_& p, const DT_ tol)
        {
          // reference hypercube: all coordinates are in range [-1,+1]
          for(int i(0); i < dim_; ++i)
          {
            // coord outside range [-1-tol, +1+tol] ?
            if((p[i] < -DT_(1)-tol) || (p[i] > DT_(1)+tol))
              return false;
          }
          return true;
        }
      };
    } // namespace Intern
  } // namespace Trafo
} // namespace FEAT