TangentOperatorBlocksView.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <MGIS/Behaviour/Variable.hxx>
#include <boost/mp11.hpp>
#include <range/v3/view/enumerate.hpp>
#include <range/v3/view/filter.hpp>
#include <range/v3/view/map.hpp>
#include <range/v3/view/zip.hpp>
 
#include "BaseLib/cpp23.h"
#include "MaterialLib/MPL/Utils/Tensor.h"
#include "MathLib/KelvinVector.h"
 
namespace MaterialLib::Solids::MFront
{
namespace detail
{
static constexpr std::array<double, 3 * 3 * 3 * 3>
    OGSMFrontTangentOperatorBlocksViewZeroes = {};
}  // namespace detail

struct OGSMFrontTangentOperatorData
{
    std::vector<double> data;
};

template <typename Gradients, typename TDynForces, typename ExtStateVars>
struct ForcesGradsCombinations
{
    static_assert(boost::mp11::mp_is_set<Gradients>::value);
    static_assert(boost::mp11::mp_is_set<TDynForces>::value);
    static_assert(boost::mp11::mp_is_set<ExtStateVars>::value);
 
    using GradientsAndExtStateVars =
        boost::mp11::mp_append<Gradients, ExtStateVars>;
 
    static_assert(boost::mp11::mp_is_set<GradientsAndExtStateVars>::value);
 
    using type = boost::mp11::
        mp_product<boost::mp11::mp_list, TDynForces, GradientsAndExtStateVars>;
};

template <int DisplacementDim, typename ForcesGradsCombinations>
class OGSMFrontTangentOperatorBlocksView
{
    static_assert(boost::mp11::mp_is_set<ForcesGradsCombinations>::value);

    static constexpr std::size_t invalid_offset_ = -1;
 
public:
    explicit OGSMFrontTangentOperatorBlocksView(
        std::vector<std::pair<mgis::behaviour::Variable,
                              mgis::behaviour::Variable>> const& to_blocks)
    {
        offsets_.fill(invalid_offset_);
 
        std::vector<bool> used_blocks(
            to_blocks.size(), false);  // checked after creation of all offsets.
 
        boost::mp11::mp_for_each<ForcesGradsCombinations>(
            [&to_blocks,
             &used_blocks,
             this]<typename Force, typename GradOrExtStateVar>(
                boost::mp11::mp_list<Force, GradOrExtStateVar>)
            {
                std::size_t data_offset = 0;
                for (auto [block, is_used] :
                     ranges::views::zip(to_blocks, used_blocks))
                {
                    auto const& [force, grad] = block;
                    if (force.name == Force::name &&
                        grad.name == GradOrExtStateVar::name)
                    {
                        auto constexpr block_idx =
                            blockIndex<Force, GradOrExtStateVar>();
                        offsets_[block_idx] = data_offset;
                        is_used = true;
                        return;
                    }
 
                    data_offset += size(force.type) * size(grad.type);
                }
            });
 
        // Indices of unused blocks.
        auto indices = ranges::views::enumerate(used_blocks) |
                       ranges::views::filter([](auto const& pair)
                                             { return !pair.second; }) |
                       ranges::views::keys;
 
        if (!indices.empty())
        {
            ERR("There are unused tangent operator blocks provided by MFront. "
                "Following blocks are unused:");
 
            for (auto const i : indices)
            {
                auto const& [force, grad] = to_blocks[i];
                ERR("\t{}/{}", force.name, grad.name);
            }
            OGS_FATAL("All tangent operator blocks must be used.");
        }
    }

    template <typename Force, typename GradOrExtStateVar>
    auto block(Force,
               GradOrExtStateVar,
               OGSMFrontTangentOperatorData const& data) const
    {
        static_assert(
            boost::mp11::mp_contains<
                ForcesGradsCombinations,
                boost::mp11::mp_list<Force, GradOrExtStateVar>>::value,
            "Requested tangent block was not created in the "
            "OGSMFrontTangentOperatorBlocksView.");
 
        constexpr auto index = blockIndex<Force, GradOrExtStateVar>();
        const auto offset = offsets_[index];
 
        constexpr auto force_size = Force::template size<DisplacementDim>();
        constexpr auto grad_size =
            GradOrExtStateVar::template size<DisplacementDim>();
 
        if constexpr (grad_size == 1 && force_size == 1)
        {
            if (offset == invalid_offset_)
            {
                return 0.0;
            }
 
            return data.data[offset];
        }
        else
        {
            constexpr auto order =
                grad_size == 1 ? Eigen::ColMajor : Eigen::RowMajor;
 
            using Result = Eigen::Map<
                const Eigen::Matrix<double, force_size, grad_size, order>>;
 
            if (offset == invalid_offset_)
            {
                return Result{
                    detail::OGSMFrontTangentOperatorBlocksViewZeroes.data()};
            }
 
            return Result{data.data.data() + offset};
        }
    }
 
private:
    static constexpr std::size_t size(mgis::behaviour::Variable::Type vt)
    {
        using VT = mgis::behaviour::Variable::Type;
 
        switch (vt)
        {
            case VT::SCALAR:
                return 1;
            case VT::STENSOR:
                return MathLib::KelvinVector::kelvin_vector_dimensions(
                    DisplacementDim);
            case VT::VECTOR:
                return DisplacementDim;
            case VT::TENSOR:
                return MaterialPropertyLib::tensorSize(DisplacementDim);
            case VT::VECTOR_1D:
            case VT::VECTOR_2D:
            case VT::VECTOR_3D:
            case VT::STENSOR_1D:
            case VT::STENSOR_2D:
            case VT::STENSOR_3D:
            case VT::TENSOR_1D:
            case VT::TENSOR_2D:
            case VT::TENSOR_3D:
            case VT::HIGHER_ORDER_TENSOR:
            case VT::ARRAY:
                break;  // Unsupported variable types
        }
 
        OGS_FATAL("Unsupported MGIS variable type {}.",
                  BaseLib::to_underlying(vt));
    }

    template <typename Force, typename GradOrExtStateVar>
    static constexpr std::size_t blockIndex()
    {
        return boost::mp11::mp_find<
            ForcesGradsCombinations,
            boost::mp11::mp_list<Force, GradOrExtStateVar>>::value;
    }

    std::array<std::size_t,
               boost::mp11::mp_size<ForcesGradsCombinations>::value>
        offsets_;
};
}  // namespace MaterialLib::Solids::MFront