LocalAssemblerTraits.h

/**
 * \file
 * \copyright
 * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#pragma once

#include "NumLib/Fem/ShapeMatrixPolicy.h"

#include <cassert>
#include <type_traits>

namespace ProcessLib
{

namespace detail
{

/*! This class makes it possible to handle Eigen's block expressions transparently,
 * both for fixed-size and dynamically allocated Eigen matrices.
 *
 * \tparam ShpPol   the shape matrix policy used
 * \tparam NNodes   the number of nodes for the FEM element
 * \tparam NodalDOF the number of d.o.f. per node.
 *                  If zero, this struct's methods will determine the size of
 *                  the generated Eigen::block's at runtime. If non-zero, the
 *                  size will be determined at compile-time.
 * \tparam Dim      global spatial dimension
 */
template<typename ShpPol, unsigned NNodes, unsigned NodalDOF, unsigned Dim>
struct LocalAssemblerTraitsFixed
{
private:
    template<int N, int M>
    using Matrix = typename ShpPol::template MatrixType<N, M>;
    template<int N>
    using Vector = typename ShpPol::template VectorType<N>;

public:
    using ShapeMatrices = typename ShpPol::ShapeMatrices;

    //! Square matrix of the given space dimensionality
    using MatrixDimDim = Matrix<Dim, Dim>;

    // TODO That only works if all process variables are single component
    //      and use the same shape functions.
    //! Local matrix for the given number of d.o.f.\ per node and number of
    //! integration points
    using LocalMatrix = Matrix<NNodes*NodalDOF, NNodes*NodalDOF>;
    //! Local vector for the given number of d.o.f.\ per node and number of
    //! integration points
    using LocalVector = Vector<NNodes*NodalDOF>;

    //! Local vector for one component of one process variable.
    //! The size is the number of nodes in the element.
    using Vector1Comp = Vector<NNodes>;

    //! Laplace matrix for the given space dimensionality and number of d.o.f
    //! per node.
    using LaplaceMatrix = Matrix<Dim*NodalDOF, Dim*NodalDOF>;


    //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<const Mat>().template block<Dim, Dim>(0u, 0u))
    >::type
    blockDimDim(Mat const& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows==Dim && ncols==Dim);
        (void) nrows; (void) ncols;
        return mat.template block<Dim, Dim>(top, left);
    }
    //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
    >::type
    blockDimDim(Mat const& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows == ncols);
        return mat.block(top, left, nrows, ncols);
    }
    //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<Mat>().template block<Dim, Dim>(0u, 0u))
    >::type
    blockDimDim(Mat& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows==Dim && ncols==Dim);
        (void) nrows; (void) ncols;
        return mat.template block<Dim, Dim>(top, left);
    }
    //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
    >::type
    blockDimDim(Mat& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows == ncols);
        return mat.block(top, left, nrows, ncols);
    }

    //! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<const Mat>().template block<NNodes, NNodes>(0u, 0u))
    >::type
    blockShpShp(Mat const& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows==NNodes && ncols==NNodes);
        (void) nrows; (void) ncols;
        return mat.template block<NNodes, NNodes>(top, left);
    }
    //! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
    >::type
    blockShpShp(Mat const& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows == ncols);
        return mat.block(top, left, nrows, ncols);
    }
    //! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<Mat>().template block<NNodes, NNodes>(0u, 0u))
    >::type
    blockShpShp(Mat& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows==NNodes && ncols==NNodes);
        (void) nrows; (void) ncols;
        return mat.template block<NNodes, NNodes>(top, left);
    }
    //! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
    template<typename Mat>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
    >::type
    blockShpShp(Mat& mat,
                unsigned top, unsigned left, unsigned nrows, unsigned ncols)
    {
        assert(nrows == ncols);
        return mat.block(top, left, nrows, ncols);
    }

    //! Get a block \c NNodes x 1 starting at the \c top'th row.
    template<typename Vec>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<const Vec>().template block<NNodes, 1>(0u, 0u))
    >::type
    blockShp(Vec const& vec, unsigned top, unsigned nrows)
    {
        assert(nrows==NNodes);
        (void) nrows;
        return vec.template block<NNodes, 1>(top, 0);
    }
    //! Get a block \c NNodes x 1 starting at the \c top'th row.
    template<typename Vec>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<const Vec>().block(0u, 0u, 0u, 0u))
    >::type
    blockShp(Vec const& vec, unsigned top, unsigned nrows)
    {
        return vec.block(top, 0, nrows, 1);
    }
    //! Get a block \c NNodes x 1 starting at the \c top'th row.
    template<typename Vec>
    static
    typename std::enable_if<NodalDOF != 0,
        decltype(std::declval<Vec>().template block<NNodes, 1>(0u, 0u))
    >::type
    blockShp(Vec& vec, unsigned top, unsigned nrows)
    {
        assert(nrows==NNodes);
        (void) nrows;
        return vec.template block<NNodes, 1>(top, 0);
    }
    //! Get a block \c NNodes x 1 starting at the \c top'th row.
    template<typename Vec>
    static
    typename std::enable_if<NodalDOF == 0,
        decltype(std::declval<Vec>().block(0u, 0u, 0u, 0u))
    >::type
    blockShp(Vec& vec, unsigned top, unsigned nrows)
    {
        return vec.block(top, 0, nrows, 1);
    }
};

} // namespace detail


#ifndef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES

template<typename ShpPol, unsigned NNodes, unsigned NodalDOF, unsigned Dim>
using LocalAssemblerTraits = detail::LocalAssemblerTraitsFixed<ShpPol, NNodes, NodalDOF, Dim>;

static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
        "Inconsistent use of the macro OGS_EIGEN_DYNAMIC_SHAPE_MATRICES."
        " Maybe you forgot to include some header file.");
#else

template<typename ShpPol, unsigned NNodes, unsigned NodalDOF, unsigned Dim>
using LocalAssemblerTraits = detail::LocalAssemblerTraitsFixed<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>;

static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
        "Inconsistent use of the macro OGS_EIGEN_DYNAMIC_SHAPE_MATRICES."
        " Maybe you forgot to include some header file.");
#endif

}  // namespace ProcessLib