corner_angle.hpp

#ifndef NIHU_CORNER_ANGLE_HPP_INCLUDED
#define NIHU_CORNER_ANGLE_HPP_INCLUDED
 
#include "element.hpp"
#include "shapeset.hpp"
 
#include "boost/math/constants/constants.hpp"
 
namespace NiHu 
{
 
 
template <class Elem, class NSet, unsigned dimension = Elem::domain_t::dimension>
class corner_angle_computer;
 
template <class Elem, class Nset>
class corner_angle_computer<Elem, Nset, 1>
{
    typedef Elem elem_t;
    typedef Nset nset_t;
public:
    static double eval(elem_t const&, unsigned nset_corner_idx)
    {
        // Position DOF can be 0 (corner) or 1 (inside)
        return 0.5 * (1 + nset_t::position_dof(nset_corner_idx));
    }
};
 
template <class Elem, class Nset>
class corner_angle_computer<Elem, Nset, 2>
{
    typedef Elem elem_t;
    typedef Nset nset_t;
 
    typedef typename elem_t::domain_t domain_t;
 
    static unsigned const num_domain_corners = domain_t::num_corners;
 
public:
    static double eval(Elem const& elem, unsigned nset_corner_idx)
    {
        auto position_dof = nset_t::position_dof(nset_corner_idx);
 
        if (position_dof == 0) // Nset corner is in domain corner
        {
            int domain_corner_idx = nset_t::node_to_domain_corner(nset_corner_idx);
 
            auto cp = domain_t::get_corner((domain_corner_idx + num_domain_corners-1) % num_domain_corners);
            auto c0 = domain_t::get_corner(domain_corner_idx);
            auto cn = domain_t::get_corner((domain_corner_idx + 1)%num_domain_corners);
 
            auto dcp = (cp - c0);
            auto dcn = (cn - c0);
 
            auto dx = elem.get_dx(nset_t::corner_at(nset_corner_idx));
 
            auto dx_p = (dx * dcp).normalized();
            auto dx_n = (dx * dcn).normalized();
 
                    
            double angle = std::acos(dx_p.dot(dx_n));
 
            return angle / (2.0 * boost::math::double_constants::pi);
            
        }
        else
        {
            // Nset corner is on edge (return 0.5), or in domain (return 1.0)
            return 0.5 * position_dof;
        }
    }
};
 
} // namespace NiHu
 
#endif /* NIHU_CORNER_ANGLE_HPP_INCLUDED */