discretizationerror.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_DISCRETIZATIONERROR_HH
#define DUNE_FUFEM_DISCRETIZATIONERROR_HH
 
#include <dune/common/dotproduct.hh>
 
#include <dune/geometry/quadraturerules.hh>
 
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
namespace Dune
{
    template<int ROWS, int COLS>
    double dot(const FieldMatrix<double,ROWS,COLS>& a, const FieldMatrix<double,ROWS,COLS>& b)
    {
        double result = 0;
        for (std::size_t i=0; i<a.rows; ++i)
            for (std::size_t j=0; j<a.cols; ++j)
                result += a[i][j] * b[i][j];
 
        return result;
    }
}
 
namespace Dune::Fufem
{
 
class DiscretizationError
{
public:
 
    template <class Function1, class Function2>
    static double computeL2DifferenceSquared(const Function1& f1, const Function2 f2, QuadratureRuleKey quadKey)
    {
        // The quantity to be computed
        double differenceSquared = 0;
 
        const auto gridView = f1.entitySet().gridView();
        const int dim = decltype(gridView)::dimension;
 
        // First function always has to be defined on a grid
        auto localFunction1 = localFunction(f1);
 
        using Element = typename decltype(gridView)::template Codim<0>::Entity;
        if constexpr (requires (const Element& e) {localFunction(f2).bind(e);} )
        {   // If second function can be evaluated in local coordinates
            auto localFunction2 = localFunction(f2);
 
            for (auto&& element : elements(gridView))
            {
                localFunction1.bind(element);
                localFunction2.bind(element);
 
                // Get quadrature formula
                quadKey.setGeometryType(element.type());
                const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
 
                for (auto&& qp : quad)
                {
                    const auto integrationElement = element.geometry().integrationElement(qp.position());
 
                    // Evaluate functions
                    const auto value1 = localFunction1(qp.position());
                    const auto value2 = localFunction2(qp.position());
 
                    // integrate squared difference
                    differenceSquared += Dune::dot(value1 - value2, value1 - value2) * qp.weight() * integrationElement;
                }
            }
        }
        else
        {   // Second function has to be evaluated in global coordinates
            for (auto&& element : elements(gridView))
            {
                localFunction1.bind(element);
 
                // Get quadrature formula
                quadKey.setGeometryType(element.type());
                const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
 
                for (auto&& qp : quad)
                {
                    const auto integrationElement = element.geometry().integrationElement(qp.position());
 
                    // Evaluate functions
                    const auto value1 = localFunction1(qp.position());
                    const auto value2 = f2(element.geometry().global(qp.position()));
 
                    // integrate squared difference
                    differenceSquared += Dune::dot(value1 - value2, value1 - value2) * qp.weight() * integrationElement;
                }
            }
        }
 
        return differenceSquared;
    }
};
 
}  // namespace Dune::Fufem
 
#endif