inverse_mapping.hpp

Go to the documentation of this file.

 
#ifndef INVERSE_MAPPING_HPP_INCLUDED
#define INVERSE_MAPPING_HPP_INCLUDED
 
#include <Eigen/Dense>
 
#include "element.hpp"
#include "shapeset.hpp"
 
#include <string>
 
namespace NiHu
{
template <class Elem>
class inverse_mapping;
 
template <class LSet, class Scalar>
class inverse_mapping<surface_element<LSet, Scalar> >
{
    typedef surface_element<LSet, Scalar> elem_t;
    typedef typename elem_t::xi_t xi_t;
    typedef typename elem_t::x_t x_t;
    typedef typename elem_t::domain_t domain_t;
    static unsigned const xi_dim = domain_t::dimension;
    static unsigned const N = xi_dim + 1;
    typedef Eigen::Matrix<Scalar, N, 1> result_t;
 
public:
    inverse_mapping(elem_t const &elem, bool constrained_in_elem)
        : m_elem(elem)
        , m_constrained_in_elem(constrained_in_elem)
    {
    }
 
    bool eval(x_t const &x0, double tol, unsigned max_iter)
    {
        // tolerance in intrinsic coordinates (only for tria and quad elems)
        // \todo xi_tol only for N = 3 surface elements
        double xi_tol = tol / std::sqrt(m_elem.get_normal(domain_t::get_center()).norm());
 
        // initial guess
        m_res.setZero();
 
        // Newton Raphson derivative matrix
        Eigen::Matrix<double, N, N> df;
 
        xi_t xi_prev = xi_t::Zero();
 
        // Newton Raphson iterations
        for (m_iter = 0; m_iter < max_iter; ++m_iter)
        {
            // intrinsic coordinates within the element
            auto xi = m_res.topRows(xi_dim);
            double zeta = m_res(xi_dim);
 
            // compute physical coordinates of guess
            auto y = m_elem.get_x(xi);
            x_t n = m_elem.get_normal(xi);
 
            double jac = n.norm();
            if (jac == 0.0)
                throw std::runtime_error("Element " + std::to_string(m_elem.get_id()(0)) + " is overdistorted, and inverse mapping cannot be performed for nearly singular integration.");
            double sqrtjac = std::sqrt(n.norm());
 
            n /= sqrtjac;
 
            auto x_guess = y + n * zeta;
 
            // check if converged
            auto f = x_guess - x0;
            m_error = f.norm();
            if (m_error < tol)
                return true;
 
            // compute derivatives
            auto dy = m_elem.get_dx(xi);
            auto ddy = m_elem.get_ddx(xi);
 
 
            if (N == 3)
            {
                auto const &dyxi = dy.col(shape_derivative_index::dXI);
                auto const &dyeta = dy.col(shape_derivative_index::dETA);
 
                auto const &ddyxixi = ddy.col(shape_derivative_index::dXIXI);
                auto const &ddyxieta = ddy.col(shape_derivative_index::dXIETA);
                auto const &ddyetaeta = ddy.col(shape_derivative_index::dETAETA);
 
                x_t dnxi = (ddyxixi.cross(dyeta) + dyxi.cross(ddyxieta)) / sqrtjac;
                x_t dneta = (ddyxieta.cross(dyeta) + dyxi.cross(ddyetaeta)) / sqrtjac;
 
                df.col(0) = dyxi + dnxi * zeta;
                df.col(1) = dyeta + dneta * zeta;
                df.col(2) = n;
            }
            else
                throw std::logic_error("Inverse mapping is unimplemented for this element type");
 
            m_res = m_res - df.colPivHouseholderQr().solve(f);
 
            if (m_constrained_in_elem)
            {
                m_res.topRows(xi_dim) = domain_t::constrain_inside(m_res.topRows(xi_dim));
                xi_t delta = m_res.topRows(xi_dim) - xi_prev;
                if (delta.norm() < xi_tol)
                    return true;
                xi_prev = m_res.topRows(xi_dim);
            }
        }
 
        return false;
    }
 
    result_t const &get_result() const
    {
        return m_res;
    }
 
    unsigned get_iter() const
    {
        return m_iter;
    }
 
    double get_error() const
    {
        return m_error;
    }
 
private:
    elem_t m_elem;
    bool m_constrained_in_elem;
    result_t m_res;
    unsigned m_iter;
    double m_error;
};
 
} // end of namespace NiHu
 
#endif // INVERSE_MAPPING_HPP_INCLUDED