LIE/SmallDeformation/LocalAssembler/LocalDataInitializer.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 <functional>
#include <memory>
#include <type_traits>
#include <typeindex>
#include <typeinfo>
#include <unordered_map>
 
#include "MeshLib/Elements/Elements.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
#include "NumLib/Fem/Integration/IntegrationMethodRegistry.h"
#include "ProcessLib/Utils/EnabledElements.h"
 
namespace ProcessLib
{
namespace LIE
{
namespace SmallDeformation
{
template <typename LocalAssemblerInterface,
          template <typename, int> class LocalAssemblerDataMatrix,
          template <typename, int> class LocalAssemblerDataMatrixNearFracture,
          template <typename, int> class LocalAssemblerDataFracture,
          int DisplacementDim, typename... ConstructorArgs>
class LocalDataInitializer final
{
    struct IsElementEnabled
    {
        template <typename ElementTraits>
        constexpr bool operator()(ElementTraits*) const
        {
            if constexpr (DisplacementDim < ElementTraits::ShapeFunction::DIM)
            {
                return false;
            }
 
            // exclude 0D elements
            return ElementTraits::Element::dimension >= 1;
        }
    };
 
public:
    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
 
    explicit LocalDataInitializer(
        NumLib::LocalToGlobalIndexMap const& dof_table,
        NumLib::IntegrationOrder const integration_order)
        : _dof_table(dof_table)
    {
        // REMARKS: At the moment, only a 2D mesh with 1D elements are
        // supported.
 
        using EnabledElementTraits =
            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
                std::declval<IsElementEnabled>()));
 
        BaseLib::TMP::foreach<EnabledElementTraits>(
            [this, integration_order]<typename ET>(ET*)
            {
                using MeshElement = typename ET::Element;
                using ShapeFunction = typename ET::ShapeFunction;
 
                _builder[std::type_index(typeid(MeshElement))] =
                    makeLocalAssemblerBuilder<ShapeFunction>(integration_order);
            });
    }

    LADataIntfPtr operator()(std::size_t const id,
                             MeshLib::Element const& mesh_item,
                             ConstructorArgs&&... args) const
    {
        auto const type_idx = std::type_index(typeid(mesh_item));
        auto const it = _builder.find(type_idx);
 
        if (it == _builder.end())
        {
            OGS_FATAL(
                "You are trying to build a local assembler for an unknown mesh "
                "element type ({:s})."
                " Maybe you have disabled this mesh element type in your build "
                "configuration, or a mesh element order does not match shape "
                "function order given in the project file.",
                type_idx.name());
        }
 
        auto const n_local_dof = _dof_table.getNumberOfElementDOF(id);
        auto const n_global_components =
            _dof_table.getNumberOfElementComponents(id);
        auto const varIDs = _dof_table.getElementVariableIDs(id);
 
        std::vector<unsigned> dofIndex_to_localIndex;
        if (mesh_item.getDimension() < DisplacementDim ||
            n_global_components > DisplacementDim)
        {
            dofIndex_to_localIndex.resize(n_local_dof);
            unsigned dof_id = 0;
            unsigned local_id = 0;
            for (auto i : varIDs)
            {
                for (int j = 0; j < _dof_table.getNumberOfVariableComponents(i);
                     j++)
                {
                    auto const& ms = _dof_table.getMeshSubset(i, j);
                    auto const mesh_id = ms.getMeshID();
                    for (unsigned k = 0; k < mesh_item.getNumberOfNodes(); k++)
                    {
                        MeshLib::Location l(mesh_id,
                                            MeshLib::MeshItemType::Node,
                                            getNodeIndex(mesh_item, k));
                        auto global_index = _dof_table.getGlobalIndex(l, i, j);
                        if (global_index != NumLib::MeshComponentMap::nop)
                        {
                            dofIndex_to_localIndex[dof_id++] = local_id;
                        }
                        local_id++;
                    }
                }
            }
        }
 
        return it->second(mesh_item, varIDs.size(), n_local_dof,
                          dofIndex_to_localIndex,
                          std::forward<ConstructorArgs>(args)...);
    }
 
private:
    using LADataBuilder = std::function<LADataIntfPtr(
        MeshLib::Element const& e,
        std::size_t const n_variables,
        std::size_t const local_matrix_size,
        std::vector<unsigned> const& dofIndex_to_localIndex,
        ConstructorArgs&&...)>;

    template <typename ShapeFunction>
    static LADataBuilder makeLocalAssemblerBuilder(
        NumLib::IntegrationOrder const integration_order)
    {
        return [integration_order](
                   MeshLib::Element const& e,
                   std::size_t const n_variables,
                   std::size_t const local_matrix_size,
                   std::vector<unsigned> const& dofIndex_to_localIndex,
                   ConstructorArgs&&... args)
        {
            auto const& integration_method = NumLib::IntegrationMethodRegistry::
                template getIntegrationMethod<
                    typename ShapeFunction::MeshElement>(integration_order);
 
            if (e.getDimension() == DisplacementDim)
            {
                if (dofIndex_to_localIndex.empty())
                {
                    return LADataIntfPtr{
                        new LocalAssemblerDataMatrix<ShapeFunction,
                                                     DisplacementDim>{
                            e, local_matrix_size, integration_method,
                            std::forward<ConstructorArgs>(args)...}};
                }
 
                return LADataIntfPtr{
                    new LocalAssemblerDataMatrixNearFracture<ShapeFunction,
                                                             DisplacementDim>{
                        e, n_variables, local_matrix_size,
                        dofIndex_to_localIndex, integration_method,
                        std::forward<ConstructorArgs>(args)...}};
            }
            return LADataIntfPtr{
                new LocalAssemblerDataFracture<ShapeFunction, DisplacementDim>{
                    e, n_variables, local_matrix_size, dofIndex_to_localIndex,
                    integration_method,
                    std::forward<ConstructorArgs>(args)...}};
        };
    }

    std::unordered_map<std::type_index, LADataBuilder> _builder;
 
    NumLib::LocalToGlobalIndexMap const& _dof_table;
};
 
}  // namespace SmallDeformation
}  // namespace LIE
}  // namespace ProcessLib