TemplateIsoparametric.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 <cassert>
#include <numbers>
 
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Node.h"
#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
 
namespace NumLib
{
template <class ShapeFunctionType_, class ShapeMatrixTypes_>
class TemplateIsoparametric
{
public:
    using ShapeFunctionType = ShapeFunctionType_;

    using ShapeMatrices = typename ShapeMatrixTypes_::ShapeMatrices;

    using NaturalCoordsMappingType =
        NaturalCoordinatesMapping<ShapeFunctionType, ShapeMatrices>;

    TemplateIsoparametric() : _ele(nullptr) {}

    explicit TemplateIsoparametric(const MeshLib::Element& e) : _ele(&e) {}

    const MeshLib::Element* getMeshElement() const { return _ele; }

    void setMeshElement(const MeshLib::Element& e) { this->_ele = &e; }
    void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
                               const unsigned global_dim,
                               bool is_axially_symmetric) const
    {
        NaturalCoordsMappingType::computeShapeMatrices(*_ele, natural_pt, shape,
                                                       global_dim);
        computeIntegralMeasure(is_axially_symmetric, shape);
    }

    template <ShapeMatrixType T_SHAPE_MATRIX_TYPE>
    void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
                               const unsigned global_dim,
                               bool is_axially_symmetric) const
    {
        NaturalCoordsMappingType::template computeShapeMatrices<
            T_SHAPE_MATRIX_TYPE>(*_ele, natural_pt, shape, global_dim);
        computeIntegralMeasure(is_axially_symmetric, shape);
    }

    double interpolateZerothCoordinate(
        typename ShapeMatrices::ShapeType const& N) const
    {
        auto* nodes = _ele->getNodes();
        typename ShapeMatrices::ShapeType rs(N.size());
        for (int i = 0; i < rs.size(); ++i)
        {
            rs[i] = (*nodes[i])[0];
        }
        auto const r = N.dot(rs);
        return r;
    }

    std::array<double, 3> interpolateCoordinates(
        typename ShapeMatrices::ShapeType const& N) const
    {
        auto* nodes = _ele->getNodes();
        typename ShapeMatrices::ShapeType rs(N.size());
        std::array<double, 3> interpolated_coords;
        for (int d = 0; d < 3; ++d)
        {
            for (int i = 0; i < rs.size(); ++i)
            {
                rs[i] = (*nodes[i])[d];
            }
            interpolated_coords[d] = N.dot(rs);
        }
        return interpolated_coords;
    }
 
private:
    void computeIntegralMeasure(bool is_axially_symmetric,
                                ShapeMatrices& shape) const
    {
        if (!is_axially_symmetric)
        {
            shape.integralMeasure = 1.0;
            return;
        }
 
        // Note: If an integration point is located at the rotation axis, r will
        // be zero which might lead to problems with the assembled equation
        // system.
        // E.g., for triangle elements, if an integration point is
        // located at edge of the unit triangle, it is possible that
        // r becomes zero.
        auto const r = interpolateZerothCoordinate(shape.N);
        shape.integralMeasure = 2. * std::numbers::pi * r;
    }
 
    const MeshLib::Element* _ele;
};

template <typename ShapeFunction, typename ShapeMatricesType>
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>
createIsoparametricFiniteElement(MeshLib::Element const& e)
{
    using FemType =
        NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
 
    return FemType{e};
}
}  // namespace NumLib