lib_shape.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2019  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2019  Peter Rucz <rucz@hit.bme.hu>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
#ifndef LIB_SHAPE_HPP_INCLUDED
#define LIB_SHAPE_HPP_INCLUDED
 
#include "../core/shapeset.hpp"
#include "lib_domain.hpp"
 
#include "line_1_gauss_shape_set.hpp"
#include "quad_1_gauss_shape_set.hpp"
#include "tria_1_gauss_shape_set.hpp"
 
 
namespace NiHu
{
 
typedef constant_shape_set<line_domain> line_0_shape_set;
typedef constant_shape_set<tria_domain> tria_0_shape_set;
typedef constant_shape_set<quad_domain> quad_0_shape_set;
typedef constant_shape_set<brick_domain> brick_0_shape_set;
 
 
typedef isoparam_shape_set<line_domain> line_1_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct jacobian_order<line_1_shape_set>
    {
        enum { value = 0 };
    };
 
    template <>
    struct shape_complexity<line_1_shape_set, 0> : matrix_function_complexity::general {};
 
    template <>
    struct shape_complexity<line_1_shape_set, 1> : matrix_function_complexity::constant {};
 
    template <>
    struct shape_complexity<line_1_shape_set, 2> : matrix_function_complexity::zero {};
}
 
template<>
class shape_function<line_1_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<line_1_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<line_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            1.0 - xi[0],
            1.0 + xi[0]
            ).finished() / 2.0;
    }
};
 
template<>
class shape_function<line_1_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<line_1_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<line_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return (shape_t() <<
            -0.5,
            +0.5
            ).finished();
    }
};
 
 
template<>
class shape_function<line_1_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<line_1_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<line_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return shape_t::Zero();
    }
};
 
 
typedef isoparam_shape_set<tria_domain> tria_1_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct jacobian_order<tria_1_shape_set>
    {
        enum { value = 0 };
    };
 
    template <>
    struct shape_complexity<tria_1_shape_set, 0> : matrix_function_complexity::general {};
 
    template <>
    struct shape_complexity<tria_1_shape_set, 1> : matrix_function_complexity::constant {};
 
    template <>
    struct shape_complexity<tria_1_shape_set, 2> : matrix_function_complexity::zero {};
}
 
 
template<>
class shape_function<tria_1_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<tria_1_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<tria_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            1.0 - xi[0] - xi[1],
            xi[0],
            xi[1]
            ).finished();
    }
};
 
template<>
class shape_function<tria_1_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<tria_1_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<tria_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return (shape_t() <<
            -1.0, -1.0,
            +1.0, 0.0,
            0.0, +1.0
            ).finished();
    }
};
 
 
template<>
class shape_function<tria_1_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<tria_1_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<tria_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return shape_t::Zero();
    }
};
 
 
 
typedef isoparam_shape_set<quad_domain> quad_1_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct jacobian_order<quad_1_shape_set>
    {
        enum { value = 1 };
    };
 
    template <>
    struct shape_complexity<quad_1_shape_set, 0> : matrix_function_complexity::general {};
 
    template <>
    struct shape_complexity<quad_1_shape_set, 1> : matrix_function_complexity::general {};
 
    template <>
    struct shape_complexity<quad_1_shape_set, 2> : matrix_function_complexity::constant {};
}
 
 
template<>
class shape_function<quad_1_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<quad_1_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<quad_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            (1.0 - xi[0])*(1.0 - xi[1]),
            (1.0 + xi[0])*(1.0 - xi[1]),
            (1.0 + xi[0])*(1.0 + xi[1]),
            (1.0 - xi[0])*(1.0 + xi[1])
            ).finished() / 4.0;
    }
};
 
template<>
class shape_function<quad_1_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<quad_1_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<quad_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            -(1.0 - xi[1]), -(1.0 - xi[0]),
            +(1.0 - xi[1]), -(1.0 + xi[0]),
            +(1.0 + xi[1]), +(1.0 + xi[0]),
            -(1.0 + xi[1]), +(1.0 - xi[0])
            ).finished() / 4.0;
    }
};
 
template<>
class shape_function<quad_1_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<quad_1_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<quad_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return (shape_t() <<
            0.0, +.25, 0.0,
            0.0, -.25, 0.0,
            0.0, +.25, 0.0,
            0.0, -.25, 0.0
            ).finished();
    }
};
 
 
typedef isoparam_shape_set<brick_domain> brick_1_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct jacobian_order<brick_1_shape_set>
    {
        enum { value = 1 };
    };
 
    template <unsigned Order>
    struct shape_complexity<brick_1_shape_set, Order> : matrix_function_complexity::general {};
}
 
 
template<>
class shape_function<brick_1_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<brick_1_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<brick_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            (1.0 - xi[0])*(1.0 - xi[1])*(1.0 - xi[2]),
            (1.0 + xi[0])*(1.0 - xi[1])*(1.0 - xi[2]),
            (1.0 + xi[0])*(1.0 + xi[1])*(1.0 - xi[2]),
            (1.0 - xi[0])*(1.0 + xi[1])*(1.0 - xi[2]),
            (1.0 - xi[0])*(1.0 - xi[1])*(1.0 + xi[2]),
            (1.0 + xi[0])*(1.0 - xi[1])*(1.0 + xi[2]),
            (1.0 + xi[0])*(1.0 + xi[1])*(1.0 + xi[2]),
            (1.0 - xi[0])*(1.0 + xi[1])*(1.0 + xi[2])
            ).finished() / 8.0;
    }
};
 
template<>
class shape_function<brick_1_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<brick_1_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<brick_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            (-1.0)*(1.0 - xi[1])*(1.0 - xi[2]), (1.0 - xi[0])*(-1.0)*(1.0 - xi[2]), (1.0 - xi[0])*(1.0 - xi[1])*(-1.0),
            (+1.0)*(1.0 - xi[1])*(1.0 - xi[2]), (1.0 + xi[0])*(-1.0)*(1.0 - xi[2]), (1.0 + xi[0])*(1.0 - xi[1])*(-1.0),
            (+1.0)*(1.0 + xi[1])*(1.0 - xi[2]), (1.0 + xi[0])*(+1.0)*(1.0 - xi[2]), (1.0 + xi[0])*(1.0 + xi[1])*(-1.0),
            (-1.0)*(1.0 + xi[1])*(1.0 - xi[2]), (1.0 - xi[0])*(+1.0)*(1.0 - xi[2]), (1.0 - xi[0])*(1.0 + xi[1])*(-1.0),
            (-1.0)*(1.0 - xi[1])*(1.0 + xi[2]), (1.0 - xi[0])*(-1.0)*(1.0 + xi[2]), (1.0 - xi[0])*(1.0 - xi[1])*(+1.0),
            (+1.0)*(1.0 - xi[1])*(1.0 + xi[2]), (1.0 + xi[0])*(-1.0)*(1.0 + xi[2]), (1.0 + xi[0])*(1.0 - xi[1])*(+1.0),
            (+1.0)*(1.0 + xi[1])*(1.0 + xi[2]), (1.0 + xi[0])*(+1.0)*(1.0 + xi[2]), (1.0 + xi[0])*(1.0 + xi[1])*(+1.0),
            (-1.0)*(1.0 + xi[1])*(1.0 + xi[2]), (1.0 - xi[0])*(+1.0)*(1.0 + xi[2]), (1.0 - xi[0])*(1.0 + xi[1])*(+1.0)
            ).finished() / 8.0;
    }
};
 
template<>
class shape_function<brick_1_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<brick_1_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<brick_1_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            0.0, 1.0 - xi[2], 1.0 - xi[1], 0.0, 1.0 - xi[0], 0.0,
            0.0, xi[2] - 1.0, xi[1] - 1.0, 0.0, xi[0] + 1.0, 0.0,
            0.0, 1.0 - xi[2], -xi[1] - 1.0, 0.0, -xi[0] - 1.0, 0.0,
            0.0, xi[2] - 1.0, xi[1] + 1.0, 0.0, xi[0] - 1.0, 0.0,
            0.0, xi[2] + 1.0, xi[1] - 1.0, 0.0, xi[0] - 1.0, 0.0,
            0.0, -xi[2] - 1.0, 1.0 - xi[1], 0.0, -xi[0] - 1.0, 0.0,
            0.0, xi[2] + 1.0, xi[1] + 1.0, 0.0, xi[0] + 1.0, 0.0,
            0.0, -xi[2] - 1.0, -xi[1] - 1.0, 0.0, 1.0 - xi[0], 0.0
            ).finished() / 8.0;
    }
};
 
 
// Forward declaration
class line_2_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct domain<line_2_shape_set> : line_domain {};
 
    template <>
    struct num_nodes<line_2_shape_set>
    {
        enum { value = 3 };
    };
 
    template <>
    struct polynomial_order<line_2_shape_set>
    {
        enum { value = 2 };
    };
 
    template <>
    struct jacobian_order<line_2_shape_set>
    {
        enum { value = 1 };
    };
 
    template <unsigned Order>
    struct shape_complexity<line_2_shape_set, Order>
    {
        typedef matrix_function_complexity::general type;
    };
 
    template <>
    struct shape_complexity<line_2_shape_set, 2>
    {
        typedef matrix_function_complexity::constant type;
    };
 
    template <>
    struct position_dof_vector<line_2_shape_set>
    {
        typedef tmp::vector<dof0, dof1, dof0> type;
    };
 
    template <>
    struct corner_index_vector<line_2_shape_set>
    {
        typedef tmp::vector<corner_index<0>, not_corner, corner_index<0> > type;
    };
}
 
class line_2_shape_set : public shape_set_base<line_2_shape_set>
{
public:
    typedef shape_set_base<line_2_shape_set> base_t;
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
    static corner_iterator_t corner_begin_impl()
    {
        return m_corners.cbegin();
    }
 
protected:
    static const corners_t m_corners;
};
 
 
template<>
class shape_function<line_2_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<line_2_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<line_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            -xi[0]*(1.0 - xi[0]) / 2.0,
            1.0 - xi[0] * xi[0],
            xi[0] * (1.0 + xi[0]) / 2.0
            ).finished();
    }
};
 
template<>
class shape_function<line_2_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<line_2_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<line_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return (shape_t() <<
            xi[0] - 0.5,
            -2.0*xi[0],
            xi[0] + 0.5
            ).finished();
    }
};
 
template<>
class shape_function<line_2_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<line_2_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<line_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return shape_t(1.0, -2.0, 1.0);
    }
};
 
 
// Forward declaration
class tria_2_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct domain<tria_2_shape_set> : tria_domain {};
 
    template <>
    struct num_nodes<tria_2_shape_set>
    {
        enum { value = 6 };
    };
 
    template <>
    struct polynomial_order<tria_2_shape_set>
    {
        enum { value = 2 };
    };
 
    template <>
    struct jacobian_order<tria_2_shape_set>
    {
        enum { value = 2 };
    };
 
    template <unsigned Order>
    struct shape_complexity<tria_2_shape_set, Order>
    {
        typedef matrix_function_complexity::general type;
    };
 
    template <>
    struct position_dof_vector<tria_2_shape_set>
    {
        typedef tmp::vector<dof0, dof1, dof0, dof1, dof0, dof1> type;
    };
 
    template <>
    struct corner_index_vector<tria_2_shape_set>
    {
        typedef tmp::vector<corner_index<0>, not_corner, corner_index<1>, not_corner, corner_index<2>, not_corner> type;
    };
}
 
 
class tria_2_shape_set : public shape_set_base<tria_2_shape_set>
{
public:
    typedef shape_set_base<tria_2_shape_set> base_t;
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
 
    static corner_iterator_t corner_begin_impl()
    {
        return m_corners.cbegin();
    }
 
protected:
    static const corners_t m_corners;
};
 
template<>
class shape_function<tria_2_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<tria_2_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<tria_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1];
        return (shape_t() <<
            (eta + xi - 1.0)*(2.0*eta + 2.0*xi - 1.0),
            -4.0*xi*(eta + xi - 1),
            xi*(2.0*xi - 1.0),
            4.0*eta*xi,
            eta*(2.0*eta - 1.0),
            -4.0*eta*(eta + xi - 1.0)
            ).finished();
    }
};
 
template<>
class shape_function<tria_2_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<tria_2_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<tria_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1];
        return (shape_t() <<
            4.0*eta + 4 * xi - 3, 4.0*eta + 4 * xi - 3.0,
            4.0 - 8 * xi - 4 * eta, -4.0*xi,
            4.0*xi - 1.0, 0.0,
            4.0*eta, 4.0*xi,
            0.0, 4.0*eta - 1.0,
            -4.0*eta, 4.0 - 4.0*xi - 8.0*eta
            ).finished();
    }
};
 
template<>
class shape_function<tria_2_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<tria_2_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<tria_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return (shape_t() <<
            4.0, 4.0, 4.0,
            -8.0, -4.0, 0.0,
            4.0, 0.0, 0.0,
            0.0, 4.0, 0.0,
            0.0, 0.0, 4.0,
            0.0, -4.0, -8.0
            ).finished();
    }
};
 
 
// Forward declaration
class quad_2_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct domain<quad_2_shape_set> : quad_domain {};
 
    template <>
    struct num_nodes<quad_2_shape_set>
    {
        enum { value = 9 };
    };
 
    template <>
    struct polynomial_order<quad_2_shape_set>
    {
        enum { value = 2 };
    };
 
    template <>
    struct jacobian_order<quad_2_shape_set>
    {
        enum { value = 3 };
    };
 
    template <unsigned Order>
    struct shape_complexity<quad_2_shape_set, Order>
    {
        typedef matrix_function_complexity::general type;
    };
 
    template <>
    struct position_dof_vector<quad_2_shape_set>
    {
        typedef tmp::vector<dof0, dof1, dof0, dof1, dof0, dof1, dof0, dof1, dof2> type;
    };
 
    template <>
    struct corner_index_vector<quad_2_shape_set>
    {
        typedef tmp::vector<corner_index<0>, not_corner, corner_index<1>, not_corner, corner_index<2>, not_corner, corner_index<3>, not_corner, not_corner> type;
    };
}
 
 
class quad_2_shape_set : public shape_set_base<quad_2_shape_set>
{
public:
    typedef shape_set_base<quad_2_shape_set> base_t;
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
    static corner_iterator_t corner_begin_impl()
    {
        return m_corners.cbegin();
    }
 
protected:
    static corners_t const m_corners;
};
 
 
template<>
class shape_function<quad_2_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<quad_2_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<quad_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1];
        auto _1mxi = 1 - xi, _1pxi = 1 + xi;
        auto _1meta = 1 - eta, _1peta = 1 + eta;
        return (shape_t() <<
            _1mxi*xi * _1meta*eta / 4.0,
            _1mxi*_1pxi * _1meta*(-eta) / 2.0,
            _1pxi*xi * _1meta*(-eta) / 4.0,
            _1pxi*xi * _1meta*_1peta / 2.0,
            _1pxi*xi * _1peta*eta / 4.0,
            _1mxi*_1pxi * _1peta*eta / 2.0,
            _1mxi*(-xi) * _1peta*eta / 4.0,
            _1mxi*(-xi) * _1meta*_1peta / 2.0,
            _1mxi*_1pxi * _1meta*_1peta
            ).finished();
    }
};
 
template<>
class shape_function<quad_2_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<quad_2_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<quad_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1], xi2 = xi*xi, eta2 = eta*eta;
        return (shape_t() <<
            eta*(2.0*xi - 1.0)*(eta - 1.0) / 4.0, xi*(2.0*eta - 1.0)*(xi - 1.0) / 4.0,
            -xi*eta*(eta - 1.0), -(xi2 - 1.0)*(2.0*eta - 1.0) / 2.0,
            eta*(2.0*xi + 1.0)*(eta - 1.0) / 4.0, xi*(2.0*eta - 1.0)*(xi + 1.0) / 4.0,
            -(2.0*xi + 1.0)*(eta2 - 1.0) / 2.0, -xi*eta*(xi + 1.0),
            eta*(2.0*xi + 1.0)*(eta + 1.0) / 4.0, xi*(2.0*eta + 1.0)*(xi + 1.0) / 4.0,
            -xi*eta*(eta + 1.0), -(xi2 - 1.0)*(2.0*eta + 1.0) / 2.0,
            eta*(2.0*xi - 1.0)*(eta + 1.0) / 4.0, xi*(2.0*eta + 1.0)*(xi - 1.0) / 4.0,
            -(2.0*xi - 1.0)*(eta2 - 1.0) / 2.0, -xi*eta*(xi - 1.0),
            2.0*xi*(eta2 - 1.0), 2.0*eta*(xi2 - 1.0)
            ).finished();
    }
};
 
template<>
class shape_function<quad_2_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<quad_2_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<quad_2_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1], xi2 = xi*xi, eta2 = eta*eta;
        return (shape_t() <<
            eta2 / 2.0 - eta / 2.0, eta*xi - xi / 2.0 - eta / 2.0 + 1.0 / 4.0, xi2 / 2.0 - xi / 2.0,
            -eta2 + eta, xi - 2 * eta*xi, 1.0 - xi2,
            eta2 / 2.0 - eta / 2.0, eta / 2.0 - xi / 2.0 + eta*xi - 1.0 / 4.0, xi2 / 2.0 + xi / 2.0,
            1.0 - eta2, -eta - 2 * eta*xi, -xi2 - xi,
            eta2 / 2.0 + eta / 2.0, eta / 2.0 + xi / 2.0 + eta*xi + 1.0 / 4.0, xi2 / 2.0 + xi / 2.0,
            -eta2 - eta, -xi - 2 * eta*xi, 1.0 - xi2,
            eta2 / 2.0 + eta / 2.0, xi / 2.0 - eta / 2.0 + eta*xi - 1.0 / 4.0, xi2 / 2.0 - xi / 2.0,
            1.0 - eta2, eta - 2 * eta*xi, -xi2 + xi,
            2 * eta2 - 2.0, 4 * eta*xi, 2 * xi2 - 2.0
            ).finished();
    }
};
 
 
// Forward declaration
class quad_28_shape_set;
 
namespace shape_set_traits
{
    template <>
    struct domain<quad_28_shape_set> : quad_domain {};
 
    template <>
    struct num_nodes<quad_28_shape_set>
    {
        enum { value = 8 };
    };
 
    template <>
    struct polynomial_order<quad_28_shape_set>
    {
        enum { value = 2 };
    };
 
    template <>
    struct jacobian_order<quad_28_shape_set>
    {
        enum { value = 3 };
    };
 
    template <unsigned Order>
    struct shape_complexity<quad_28_shape_set, Order>
    {
        typedef matrix_function_complexity::general type;
    };
 
    template <>
    struct position_dof_vector<quad_28_shape_set>
    {
        typedef tmp::vector<dof0, dof1, dof0, dof1, dof0, dof1, dof0, dof1> type;
    };
 
    template <>
    struct corner_index_vector<quad_28_shape_set>
    {
        typedef tmp::vector<corner_index<0>, not_corner, corner_index<1>, not_corner, corner_index<2>, not_corner, corner_index<3>, not_corner> type;
    };
}
 
class quad_28_shape_set : public shape_set_base<quad_28_shape_set>
{
public:
    typedef shape_set_base<quad_28_shape_set> base_t;
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
    static corner_iterator_t corner_begin_impl()
    {
        return m_corners.cbegin();
    }
 
protected:
    static corners_t const m_corners;
};
 
 
template<>
class shape_function<quad_28_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<quad_28_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<quad_28_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1], xi2 = xi*xi, eta2 = eta*eta;
        return (shape_t() <<
            -((xi - 1.0)*(eta - 1.0)*(xi + eta + 1.0)) / 4.0,
            ((xi2 - 1.0)*(eta - 1.0)) / 2.0,
            ((xi + 1.0)*(eta - 1.0)*(eta - xi + 1.0)) / 4.0,
            -((eta2 - 1.0)*(xi + 1.0)) / 2.0,
            ((xi + 1.0)*(eta + 1.0)*(xi + eta - 1.0)) / 4.0,
            -((xi2 - 1.0)*(eta + 1.0)) / 2.0,
            ((xi - 1.0)*(eta + 1.0)*(xi - eta + 1.0)) / 4.0,
            ((eta2 - 1.0)*(xi - 1.0)) / 2.0
            ).finished();
    }
};
 
template<>
class shape_function<quad_28_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<quad_28_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<quad_28_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto x = _xi[0], y = _xi[1], x2 = x*x, y2 = y*y;
        return (shape_t() <<
            -((2.0*x + y)*(y - 1.0)) / 4.0, -((x + 2 * y)*(x - 1.0)) / 4.0,
            x*(y - 1.0), (x2 - 1.0) / 2.0,
            -((2.0*x - y)*(y - 1.0)) / 4.0, -((x - 2 * y)*(x + 1.0)) / 4.0,
            (1.0 - y2) / 2.0, -y*(x + 1.0),
            ((2.0*x + y)*(y + 1.0)) / 4.0, ((x + 2 * y)*(x + 1.0)) / 4.0,
            -x*(y + 1.0), (1.0 - x2) / 2.0,
            ((2.0*x - y)*(y + 1.0)) / 4.0, ((x - 2 * y)*(x - 1.0)) / 4.0,
            (y2 - 1.0) / 2.0, y*(x - 1.0)
            ).finished();
    }
};
 
template<>
class shape_function<quad_28_shape_set, 2>
{
    typedef shape_set_traits::shape_value_type<quad_28_shape_set, 2>::type shape_t;
    typedef shape_set_traits::domain<quad_28_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi(_xi[0]), eta(_xi[1]);
        return (shape_t() <<
            .5 - eta / 2.0, 1 / 4.0 - xi / 2.0 - eta / 2.0, .5 - xi / 2.0,
            eta - 1.0, xi, 0.0,
            .5 - eta / 2.0, eta / 2.0 - xi / 2.0 - 1 / 4.0, xi / 2.0 + .5,
            0.0, -eta, -xi - 1.0,
            eta / 2.0 + .5, eta / 2.0 + xi / 2.0 + 1 / 4.0, xi / 2.0 + .5,
            -eta - 1.0, -xi, 0.0,
            eta / 2.0 + .5, xi / 2.0 - eta / 2.0 - 1 / 4.0, .5 - xi / 2.0,
            0.0, eta, xi - 1.0
            ).finished();
    }
};
 
} // namespace NiHu
 
 
#endif // LIB_SHAPE_HPP_INCLUDED