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/base_header.hpp>
#include <kernel/util/tiny_algebra.hpp>
 
namespace FEAT
{
  namespace Trafo
  {
    namespace Standard
    {
      template<typename DT_, int helper_dim, int world_dim, int sc_, int sp_, int smx_, int snx_,
        typename std::enable_if<helper_dim == world_dim+1, bool>::type = true>
      void inverse_mapping(
        Tiny::Vector<DT_, world_dim, sc_>& coeffs,
        const Tiny::Vector<DT_, world_dim, sp_>& point,
        const Tiny::Matrix<DT_, helper_dim, world_dim, smx_, snx_>& x)
        {
          // A will contain the transformation matrix
          Tiny::Matrix<DT_, world_dim, world_dim> A(DT_(0));
          for(int i(0); i < world_dim; ++i)
          {
            for(int j(0); j < world_dim; ++j)
            {
              A(i,j) = x(j+1,i) - x(0,i);
            }
          }
 
          // This will be the inverse of the transformation matrix
          Tiny::Matrix<DT_, world_dim, world_dim> Ainv(DT_(0));
          Ainv.set_inverse(A);
 
          // This is the solution of A u = point - x[0]
          coeffs = Ainv*(point - x[0]);
        }
 
      template<typename DT_, int world_dim, int sc_, int sp_, int smx_, int snx_>
      void inverse_mapping(
        Tiny::Vector<DT_, 2, sc_>& coeffs,
        const Tiny::Vector<DT_, world_dim, sp_>& point,
        const Tiny::Matrix<DT_, 2, world_dim, smx_, snx_>& x)
      {
        static_assert( (world_dim == 2 || world_dim == 3),
        "world dim has to be 2 or 3 for complementary barycentric coordinates");
 
        auto tmp = x[1]-x[0];
        DT_ sp(Tiny::dot(point - x[0],tmp));
        DT_ nsqr(Math::sqr(tmp.norm_euclid()));
        // This is omega
        coeffs[0] = sp/nsqr;
        tmp = point - (x[0] + coeffs[0]*(x[1]-x[0]));
        // This is the distance of point to the straight line defined by x[0], x[1]
        coeffs[1] = tmp.norm_euclid();
      }
 
      template<typename DT_, int sc_, int sp_, int smx_, int snx_>
      void inverse_mapping(
        Tiny::Vector<DT_, 3, sc_>& coeffs,
        const Tiny::Vector<DT_, 3, sp_>& point,
        const Tiny::Matrix<DT_, 3, 3, smx_, snx_>& x)
        {
          static constexpr int world_dim = 3;
          static constexpr int shape_dim = world_dim-1;
 
          // A will contain the transformation matrix
          Tiny::Matrix<DT_, world_dim, world_dim> A(DT_(0));
          // Fill the all rows in the first shape_dim = world_dim-1 columns
          for(int i(0); i < world_dim; ++i)
          {
            for(int j(0); j < shape_dim; ++j)
              A(i,j) = x(j+1,i) - x(0,i);
          }
 
          // The last column is the additional direction for our augmented simplex and it is orthogonal to the rest
          Tiny::Vector<DT_, world_dim> ortho;
          Tiny::orthogonal_3x2(ortho, A);
 
          // In 3d, the 2-norm of ortho is the 2d volume of the parallelogram defined by A[0], A[1]. That makes this
          // axis badly scaled if the parallelogram is either very small or very large. So we rescale ortho to unity
          // so the resulting coefficient is just the distance
          ortho.normalize();
 
          // Set the last column in A
          for(int i(0); i < world_dim; ++i)
            A(i,world_dim-1) = ortho(i);
 
          // This will be the inverse of the transformation matrix
          Tiny::Matrix<DT_, world_dim, world_dim> Ainv(DT_(0));
          Ainv.set_inverse(A);
 
          // This is the solution of A u = point - x[0]
          coeffs = Ainv*(point - x[0]);
        }
 
    } // namespace Standard
  } // namespace Trafo
} // namespace FEAT