LocalDataInitializer.h

Go to the documentation of this file.

 
#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/FiniteElement/LowerDimShapeTable.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
 
#ifndef OGS_MAX_ELEMENT_DIM
static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
#endif
 
#ifndef OGS_MAX_ELEMENT_ORDER
static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
#endif
 
// The following macros decide which element types will be compiled, i.e.
// which element types will be available for use in simulations.
 
#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
#else
#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
#endif
 
#ifdef OGS_ENABLE_ELEMENT_CUBOID
#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
#else
#define ENABLED_ELEMENT_TYPE_CUBOID 0u
#endif
 
#ifdef OGS_ENABLE_ELEMENT_PRISM
#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
#else
#define ENABLED_ELEMENT_TYPE_PRISM 0u
#endif
 
#ifdef OGS_ENABLE_ELEMENT_PYRAMID
#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
#else
#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
#endif
 
// Dependent element types.
// Faces of tets, pyramids and prisms are triangles
#define ENABLED_ELEMENT_TYPE_TRI                                       \
    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
     (ENABLED_ELEMENT_TYPE_PRISM))
// Faces of hexes, pyramids and prisms are quads
#define ENABLED_ELEMENT_TYPE_QUAD                                     \
    ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
     (ENABLED_ELEMENT_TYPE_PRISM))
 
// All enabled element types
#define OGS_ENABLED_ELEMENTS                                          \
    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
     (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
 
// Include only what is needed (Well, the conditions are not sharp).
#if OGS_ENABLED_ELEMENTS != 0
#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
#endif
 
namespace ProcessLib
{
namespace LIE
{
namespace HydroMechanics
{
template <typename LocalAssemblerInterface,
          template <typename, typename, typename, int>
          class LocalAssemblerDataMatrix,
          template <typename, typename, typename, int>
          class LocalAssemblerDataMatrixNearFracture,
          template <typename, typename, typename, int>
          class LocalAssemblerDataFracture,
          int GlobalDim, typename... ConstructorArgs>
class LocalDataInitializer final
{
public:
    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
 
    LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
                         const unsigned shapefunction_order)
        : _dof_table(dof_table)
    {
        if (shapefunction_order != 2)
        {
            OGS_FATAL(
                "The given shape function order {:d} is not supported.\nOnly "
                "shape functions of order 2 are supported.",
                shapefunction_order);
        }
            // /// Quads and Hexahedra ///////////////////////////////////
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Quad8))] =
            makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
        _builder[std::type_index(typeid(MeshLib::Quad9))] =
            makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Hex20))] =
            makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
#endif
 
        // /// Simplices ////////////////////////////////////////////////
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Tri6))] =
            makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
#endif
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Tet10))] =
            makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
#endif
 
        // /// Prisms ////////////////////////////////////////////////////
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Prism15))] =
            makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
#endif
 
        // /// Pyramids //////////////////////////////////////////////////
 
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
            makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
#endif
        // /// Lines ///////////////////////////////////
 
// /// Lines ///////////////////////////////////
#if OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
        _builder[std::type_index(typeid(MeshLib::Line3))] =
            makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
#endif
    }
 
    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 varIDs = _dof_table.getElementVariableIDs(id);
        bool const isPressureDeactivated = (varIDs.front() != 0);
        std::vector<int> involved_varIDs;  // including deactivated elements
        involved_varIDs.reserve(varIDs.size() + 1);
        if (isPressureDeactivated)
        {
            involved_varIDs.push_back(0);  // always pressure come in
        }
        involved_varIDs.insert(involved_varIDs.end(), varIDs.begin(),
                               varIDs.end());
 
        std::vector<unsigned> dofIndex_to_localIndex;
 
        // matrix and fracture assemblers with enrichments
        dofIndex_to_localIndex.resize(n_local_dof);
        std::vector<unsigned> vec_n_element_nodes;
        // TODO how to get the shape function order for each variable?
        vec_n_element_nodes.push_back(
            mesh_item.getNumberOfBaseNodes());  // pressure
        auto const max_varID = *std::max_element(varIDs.begin(), varIDs.end());
        for (int i = 1; i < max_varID + 1; i++)
        {
            vec_n_element_nodes.push_back(
                mesh_item.getNumberOfNodes());  // displacements
        }
 
        unsigned local_id = 0;
        unsigned dof_id = 0;
        for (unsigned i = 0; i < involved_varIDs.size(); i++)
        {
            auto const var_id = involved_varIDs[i];
            auto const n_var_comp =
                _dof_table.getNumberOfVariableComponents(var_id);
            auto const n_var_element_nodes = vec_n_element_nodes[i];
            for (int var_comp_id = 0; var_comp_id < n_var_comp; var_comp_id++)
            {
                auto const& ms = _dof_table.getMeshSubset(var_id, var_comp_id);
                auto const mesh_id = ms.getMeshID();
                for (unsigned k = 0; k < n_var_element_nodes; k++)
                {
                    MeshLib::Location l(mesh_id,
                                        MeshLib::MeshItemType::Node,
                                        getNodeIndex(mesh_item, k));
                    auto global_index =
                        _dof_table.getGlobalIndex(l, var_id, var_comp_id);
                    if (global_index != NumLib::MeshComponentMap::nop)
                    {
                        dofIndex_to_localIndex[dof_id++] = local_id;
                    }
                    local_id++;
                }
            }
        }
 
        return it->second(mesh_item, involved_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 ShapeFunctionDisplacement>
    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
        typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
 
    template <typename ShapeFunctionDisplacement,
              typename ShapeFunctionPressure>
    using LADataMatrix = LocalAssemblerDataMatrix<
        ShapeFunctionDisplacement, ShapeFunctionPressure,
        IntegrationMethod<ShapeFunctionDisplacement>, GlobalDim>;
 
    template <typename ShapeFunctionDisplacement,
              typename ShapeFunctionPressure>
    using LADataMatrixNearFracture = LocalAssemblerDataMatrixNearFracture<
        ShapeFunctionDisplacement, ShapeFunctionPressure,
        IntegrationMethod<ShapeFunctionDisplacement>, GlobalDim>;
 
    template <typename ShapeFunctionDisplacement,
              typename ShapeFunctionPressure>
    using LAFractureData = LocalAssemblerDataFracture<
        ShapeFunctionDisplacement, ShapeFunctionPressure,
        IntegrationMethod<ShapeFunctionDisplacement>, GlobalDim>;
 
    template <typename ShapeFunctionDisplacement>
    static LADataBuilder makeLocalAssemblerBuilder()
    {
        return makeLocalAssemblerBuilder<ShapeFunctionDisplacement>(
            static_cast<std::integral_constant<
                bool, (GlobalDim >= ShapeFunctionDisplacement::DIM)>*>(
                nullptr));
    }
 
    std::unordered_map<std::type_index, LADataBuilder> _builder;
 
    NumLib::LocalToGlobalIndexMap const& _dof_table;
 
private:
    template <typename ShapeFunctionDisplacement>
    static LADataBuilder makeLocalAssemblerBuilder(std::true_type* /*unused*/)
    {
        // (Lower order elements = Order(ShapeFunctionDisplacement) - 1).
        using ShapeFunctionPressure =
            typename NumLib::LowerDim<ShapeFunctionDisplacement>::type;
        return [](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) {
            if (e.getDimension() == GlobalDim)
            {
                if (n_variables == 2)
                {
                    return LADataIntfPtr{
                        new LADataMatrix<ShapeFunctionDisplacement,
                                         ShapeFunctionPressure>{
                            e, n_variables, local_matrix_size,
                            dofIndex_to_localIndex,
                            std::forward<ConstructorArgs>(args)...}};
                }
                return LADataIntfPtr{new LADataMatrixNearFracture<
                    ShapeFunctionDisplacement, ShapeFunctionPressure>{
                    e, n_variables, local_matrix_size, dofIndex_to_localIndex,
                    std::forward<ConstructorArgs>(args)...}};
            }
            return LADataIntfPtr{new LAFractureData<ShapeFunctionDisplacement,
                                                    ShapeFunctionPressure>{
                e, local_matrix_size, dofIndex_to_localIndex,
                std::forward<ConstructorArgs>(args)...}};
        };
    }
 
    template <typename ShapeFunctionDisplacement>
    static LADataBuilder makeLocalAssemblerBuilder(std::false_type* /*unused*/)
    {
        return nullptr;
    }
};
 
}  // namespace HydroMechanics
}  // namespace LIE
}  // namespace ProcessLib
 
#undef ENABLED_ELEMENT_TYPE_SIMPLEX
#undef ENABLED_ELEMENT_TYPE_CUBOID
#undef ENABLED_ELEMENT_TYPE_PYRAMID
#undef ENABLED_ELEMENT_TYPE_PRISM
#undef ENABLED_ELEMENT_TYPE_TRI
#undef ENABLED_ELEMENT_TYPE_QUAD
#undef OGS_ENABLED_ELEMENTS