Insufficient quadrature schemes

[dbeurle]

The quadrature schemes lack a degree definition, which is the important factor in deciding which rule to use for what bilinear form and weight & test function combination. For example, the mass (or mass-like) matrix integral with quadratic tetrahedrons omits a fourth order polynomial which requires a fourth order accurate scheme.

Tasks:

• [x] Add quadrature degree interface
• [x] Add high order quadrature rules
• [x] Implement methods to shape function to return poly- and mononomial degrees

Selection criteria: A necessary (but not sufficient) condition is that the order of the scheme must be k_bar + k - 2m, where k_bar is the monomial order, k is the polynomial order and m is the order of the derivative that appears in the weak form.

• [x] Implement full convergence order check for the error norm
• [x] Compare order for mass matrix
• [ ] Add rank deficiency check

After numerical experiments, an insufficient quadrature scheme for the mass matrix will yield zero, or close to zero eigenvalues which degrades the condition number of the mass matrix. This is consistent with the higher order terms in the mass matrix. Since the mass matrix can be integrated once and does not depend on internal variables, a separate scheme can be used.

To this end I propose the following changes to the numerical integration and shape function classes. Here it is important to note the construction dependency: derivative order, shape function and then an appropriate integration degree can be chosen.

The numerical_integration class will only store the coordinates, degree and the weights of the shape function and make these accessible through constant methods. It will remain templated over the number of doubles which span the space (one, two or three dimensions) and aliased to line_quadrature, surface_quadrature and volume_quadrature.

The shape_function class will be changed to be no longer coupled with the quadrature class. A type alias will be introduced as

using coordinate_type = std::tuple<int, double, double, double>;

Additional methods will be provided to the interface:

/// Single coordinate point
auto evaluate(coordinate_type const& coordinate) const -> std::pair<vector, matrix>

/// Multiple coordinate points
auto evaluate(std::vector<coordinate_type> const& coordinates) const -> std::vector<std::pair<vector, matrix>>

which evaluates the shape functions and shape function derivatives at the specified coordinates. Additionally a method that returns the local element coordinates:

auto local_coordinates() const -> std::vector<coordinate_type> const&

Alone the quadrature and the shape function classes are rather useless, so their functionality will be combined into a integral_form (name subject to change) class that handles the numerical integration. Any remaining functionality of numerical_integration will be moved into this class.

template <typename Interpolation, typename Quadrature>
class integral_form
{
public:
    /// volume, surface or line interpolation
    using interpolation_type = Interpolation;
    /// volume, surface or line quadrature
    using quadrature_type = Quadrature;

public:
    integral_form(int const derivative_order) 
        // allocate interpolation type
        // allocate quadrature using interpolation and derivative_order
        // m_values{m_interpolation->evaluate(m_quadrature.coordinates())}
    {
    }

    template <typename Expression>
    void integrate_inplace(matrix& integrand, Expression&& expression)
    {
        auto const& weights = m_quadrature->weights();

        std::size_t const size = weights.size();

        for (std::size_t index = 0; index < size; ++index)
        {
            integrand += expression(m_values[index]) * weights[index];
        }
    }

    /// if required
    auto interpolation() const -> interpolation_type const&;
    auto quadrature() const -> quadrature_type const&;

protected:
    std::unique_ptr<interpolation_type> m_interpolation;
    std::unique_ptr<quadrature_type> m_quadrature;

    std::vector<std::pair<vector, matrix>> m_values;
};