IntegrationMethodRegistry.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "IntegrationMethodRegistry.h"
 
#include <unordered_map>
 
#include "BaseLib/DemangleTypeInfo.h"
#include "BaseLib/Error.h"
#include "GaussLegendreIntegrationPolicy.h"
#include "MeshLib/Elements/Elements.h"
 
static constexpr unsigned MAX_ORDER_REGULAR = 4u;
 
template <typename MeshElement>
using IntegrationPolicy = NumLib::GaussLegendreIntegrationPolicy<MeshElement>;
 
template <typename IntegrationPolicy_>
static NumLib::GenericIntegrationMethod createGenericIntegrationMethod(
    unsigned const order)
{
    typename IntegrationPolicy_::IntegrationMethod meth{order};
    unsigned const np = meth.getNumberOfPoints();
 
    std::vector<MathLib::WeightedPoint> wps;
    wps.reserve(np);
 
    for (unsigned ip = 0; ip < np; ++ip)
    {
        wps.emplace_back(meth.getWeightedPoint(ip));
 
        if (wps.back() != meth.getWeightedPoint(ip))
        {
            throw std::runtime_error(
                "createGenericIntegrationMethod mismatch for ip=" +
                std::to_string(ip) + ", order=" + std::to_string(order) +
                ", method=" + BaseLib::typeToString(meth));
        }
    }
 
    return NumLib::GenericIntegrationMethod{order, std::move(wps)};
}
 
template <typename MeshElement>
static void putIntegrationMethodsFor(
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>&
        integration_methods_by_mesh_element_type,
    unsigned const max_order)
{
    std::vector<NumLib::GenericIntegrationMethod> integration_methods;
    integration_methods.reserve(max_order + 1);
 
    // order 0 -> no valid weighted points
    NumLib::GenericIntegrationMethod invalidMethod{0, {}};
    integration_methods.push_back(std::move(invalidMethod));
 
    for (unsigned order = 1; order <= max_order; ++order)
    {
        integration_methods.push_back(
            createGenericIntegrationMethod<IntegrationPolicy<MeshElement>>(
                order));
    }
 
    integration_methods_by_mesh_element_type.emplace(
        std::type_index(typeid(MeshElement)), std::move(integration_methods));
}
 
static void putIntegrationMethodsForDim0(
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>&
        integration_methods_by_mesh_element_type)
{
    putIntegrationMethodsFor<MeshLib::Point>(
        integration_methods_by_mesh_element_type,
        MAX_ORDER_REGULAR /* arbitrary cutoff */);
}
 
static void putIntegrationMethodsForDim1(
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>&
        integration_methods_by_mesh_element_type)
{
    putIntegrationMethodsFor<MeshLib::Line>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Line3>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
}
 
static void putIntegrationMethodsForDim2(
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>&
        integration_methods_by_mesh_element_type)
{
    putIntegrationMethodsFor<MeshLib::Quad>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Quad8>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Quad9>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Tri>(
        integration_methods_by_mesh_element_type, 4);
 
    putIntegrationMethodsFor<MeshLib::Tri6>(
        integration_methods_by_mesh_element_type, 4);
}
 
static void putIntegrationMethodsForDim3(
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>&
        integration_methods_by_mesh_element_type)
{
    putIntegrationMethodsFor<MeshLib::Hex>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Hex20>(
        integration_methods_by_mesh_element_type, MAX_ORDER_REGULAR);
 
    putIntegrationMethodsFor<MeshLib::Tet>(
        integration_methods_by_mesh_element_type, 4);
 
    putIntegrationMethodsFor<MeshLib::Tet10>(
        integration_methods_by_mesh_element_type, 4);
 
    putIntegrationMethodsFor<MeshLib::Prism>(
        integration_methods_by_mesh_element_type, 4);
 
    putIntegrationMethodsFor<MeshLib::Prism15>(
        integration_methods_by_mesh_element_type, 4);
 
    /* Note: Currently (July 22) OGS has pyramid integration schemes only up to
     * order 3 (see MathLib/Integration/GaussLegendrePyramid.h), and for order 4
     * the third order is reused (see
     * NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h).
     * I.e., a user can request 4th order integration on pyramids and OGS will
     * just run fine, albeit with a lower integration order.
     * To keep that behavior, in particular not crash if a user requests 4th
     * order pyramid integration, we set the limit to 4, here.
     */
 
    putIntegrationMethodsFor<MeshLib::Pyramid>(
        integration_methods_by_mesh_element_type, 4 /* see note above */);
 
    putIntegrationMethodsFor<MeshLib::Pyramid13>(
        integration_methods_by_mesh_element_type, 4 /* see note above */);
}
 
static std::unordered_map<std::type_index,
                          std::vector<NumLib::GenericIntegrationMethod>>
initIntegrationMethods()
{
    std::unordered_map<std::type_index,
                       std::vector<NumLib::GenericIntegrationMethod>>
        integration_methods_by_mesh_element_type;
 
    putIntegrationMethodsForDim0(integration_methods_by_mesh_element_type);
    putIntegrationMethodsForDim1(integration_methods_by_mesh_element_type);
    putIntegrationMethodsForDim2(integration_methods_by_mesh_element_type);
    putIntegrationMethodsForDim3(integration_methods_by_mesh_element_type);
 
    return integration_methods_by_mesh_element_type;
}
 
namespace NumLib::IntegrationMethodRegistry
{
GenericIntegrationMethod const& getIntegrationMethod(
    std::type_index const mesh_element_type, IntegrationOrder const order)
{
    if (order.order == 0)
    {
        // For 0D elements (points) the order does not matter, still we don't
        // want order zero there. For all other element types integration order
        // does matter.
        OGS_FATAL("An integration order of 0 is not supported.");
    }
 
    // The actual integration method registry.
    //
    // A mapping \code [mesh element type] -> [list of integration
    // methods]\endcode, where each entry in the list corresponds to integration
    // order 0, 1, 2, ...
    //
    // Having this as a static __local__ variable circumvents the static
    // initialization order fiasco we would have if this was a static __global__
    // variable. The fiasco arises, because this registry uses static global
    // data from MathLib. Currently (Feb 2022) we cannot circumvent the problem
    // with constinit due to lack of compiler support.
    static const std::unordered_map<
        std::type_index,
        std::vector<NumLib::GenericIntegrationMethod>>
        integration_methods_by_mesh_element_type = initIntegrationMethods();
 
    if (auto it =
            integration_methods_by_mesh_element_type.find(mesh_element_type);
        it != integration_methods_by_mesh_element_type.end())
    {
        auto& integration_methods = it->second;
 
        if (order.order >= integration_methods.size())
        {
            OGS_FATAL(
                "Integration order {} is not supported for mesh elements of "
                "type {}",
                order.order,
                mesh_element_type.name());
        }
 
        return integration_methods[order.order];
    }
 
    OGS_FATAL(
        "No integration methods are available for mesh elements of type {}",
        mesh_element_type.name());
}
}  // namespace NumLib::IntegrationMethodRegistry