MeshQuadrature.h

Go to the documentation of this file.

 
#ifndef PBAT_FEM_MESHQUADRATURE_H
#define PBAT_FEM_MESHQUADRATURE_H
 
#include "Concepts.h"
#include "Jacobian.h"
#include "ShapeFunctions.h"
#include "pbat/Aliases.h"
#include "pbat/common/Concepts.h"
 
#include <exception>
#include <fmt/format.h>
 
namespace pbat::fem {

template <CElement TElement, int QuadratureOrder, class TDerivedDetJe>
void ToMeshQuadratureWeights(Eigen::MatrixBase<TDerivedDetJe>& detJeThenWg)
{
    PBAT_PROFILE_NAMED_SCOPE("pbat.fem.ToMeshQuadratureWeights");
    using ScalarType     = typename TDerivedDetJe::Scalar;
    using QuadratureType = typename TElement::template QuadratureType<QuadratureOrder, ScalarType>;
    if (detJeThenWg.rows() != QuadratureType::kPoints)
    {
        throw std::invalid_argument(
            fmt::format(
                "detJeThenWg must have {} rows, but has {}",
                QuadratureType::kPoints,
                detJeThenWg.rows()));
    }
    auto const wg = common::ToEigen(QuadratureType::weights);
    detJeThenWg.array().colwise() *= wg.array();
}

template <CElement TElement, int QuadratureOrder, class TDerivedE, class TDerivedX>
auto MeshQuadratureWeights(
    Eigen::MatrixBase<TDerivedE> const& E,
    Eigen::MatrixBase<TDerivedX> const& X)
    -> Eigen::Matrix<
        typename TDerivedX::Scalar,
        TElement::template QuadratureType<QuadratureOrder, typename TDerivedX::Scalar>::kPoints,
        Eigen::Dynamic>
{
    PBAT_PROFILE_NAMED_SCOPE("pbat.fem.MeshQuadratureWeights");
    using ElementType = TElement;
    auto detJeThenWg = DeterminantOfJacobian<ElementType, QuadratureOrder>(E, X);
    ToMeshQuadratureWeights<ElementType, QuadratureOrder>(detJeThenWg);
    return detJeThenWg;
}

template <int QuadratureOrder, CMesh TMesh>
auto MeshQuadratureWeights(TMesh const& mesh) -> Eigen::Matrix<
    typename TMesh::ScalarType,
    TMesh::ElementType::template QuadratureType<QuadratureOrder, typename TMesh::ScalarType>::
        kPoints,
    Eigen::Dynamic>
{
    return MeshQuadratureWeights<typename TMesh::ElementType, QuadratureOrder>(mesh.E, mesh.X);
}

template <common::CIndex TIndex>
auto MeshQuadratureElements(TIndex nElements, TIndex nQuadPtsPerElement)
{
    using IndexVectorType = Eigen::Vector<TIndex, Eigen::Dynamic>;
    return IndexVectorType::LinSpaced(nElements, TIndex(0), nElements - 1)
        .replicate(TIndex(1), nQuadPtsPerElement)
        .transpose();
}

template <common::CIndex TIndex, class TDerivedwg>
auto MeshQuadratureElements(TIndex nElements, Eigen::DenseBase<TDerivedwg> const& wg)
{
    bool const bAreDimensionsCorrect = (wg.size() % nElements) == 0;
    if (not bAreDimensionsCorrect)
    {
        throw std::invalid_argument(
            fmt::format(
                "nElements must match the number of columns in wg, or wg's flattened size must be "
                "a multiple of nElements but got {} and {}",
                nElements,
                wg.size()));
    }
    auto const nQuadPtsPerElement = wg.size() / nElements;
    return MeshQuadratureElements(nElements, nQuadPtsPerElement);
}

template <class TDerivedE, class TDerivedwg>
auto MeshQuadratureElements(
    Eigen::DenseBase<TDerivedE> const& E,
    Eigen::DenseBase<TDerivedwg> const& wg)
{
    return MeshQuadratureElements(E.cols(), wg);
}

template <
    CElement TElement,
    int QuadratureOrder,
    common::CArithmetic TScalar = Scalar,
    common::CIndex TIndex       = Index>
auto MeshReferenceQuadraturePoints(TIndex nElements)
{
    using QuadratureType = typename TElement::template QuadratureType<QuadratureOrder, TScalar>;
    auto const Xi =
        common::ToEigen(QuadratureType::points).template bottomRows<QuadratureType::kDims>();
    return Xi.replicate(1, nElements);
}

template <int QuadratureOrder, CMesh TMesh>
auto MeshReferenceQuadraturePoints(TMesh const& mesh)
{
    using ElementType = typename TMesh::ElementType;
    using ScalarType  = typename TMesh::ScalarType;
    return MeshReferenceQuadraturePoints<ElementType, QuadratureOrder, ScalarType>(mesh.E.cols());
}
 
} // namespace pbat::fem
 
#endif // PBAT_FEM_MESHQUADRATURE_H