LagrangeSpace.h

// Class LagrangeSpace
//    This is the LagrangeSpace class. It is a class representing a Lagrange space
//    for finite element methods on a structured grid.
 
#ifndef IPPL_LAGRANGESPACE_H
#define IPPL_LAGRANGESPACE_H
 
#include <cmath>
 
#include "FEM/FEMQuadratureData.h"
#include "FEM/FiniteElementSpace.h"
 
constexpr unsigned getLagrangeNumElementDOFs(unsigned Dim, unsigned Order) {
    // needs to be constexpr pow function to work at compile time. Kokkos::pow doesn't work.
    return static_cast<unsigned>(power(static_cast<int>(Order + 1), static_cast<int>(Dim)));
}
 
namespace ippl {
 
    template <typename T, unsigned Dim, unsigned Order, typename ElementType,
              typename QuadratureType, typename FieldLHS, typename FieldRHS>
    // requires IsQuadrature<QuadratureType>
    class LagrangeSpace
        : public FiniteElementSpace<T, Dim, getLagrangeNumElementDOFs(Dim, Order), ElementType,
                                    QuadratureType, FieldLHS, FieldRHS> {
    public:
        // The number of degrees of freedom per element
        static constexpr unsigned numElementDOFs = getLagrangeNumElementDOFs(Dim, Order);
 
        // The dimension of the mesh
        static constexpr unsigned dim = FiniteElementSpace<T, Dim, numElementDOFs, ElementType,
                                                           QuadratureType, FieldLHS, FieldRHS>::dim;
 
        // The order of the Lagrange space
        static constexpr unsigned order = Order;
 
        // The number of mesh vertices per element
        static constexpr unsigned numElementVertices =
            FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FieldLHS,
                               FieldRHS>::numElementVertices;
 
        // A vector with the position of the element in the mesh in each dimension
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FieldLHS, FieldRHS>::indices_t indices_t;
 
        // A point in the global coordinate system
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FieldLHS, FieldRHS>::point_t point_t;
 
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FieldLHS, FieldRHS>::vertex_points_t vertex_points_t;
 
        // Field layout type for domain decomposition info
        typedef FieldLayout<Dim> Layout_t;
 
        // View types
        typedef typename detail::ViewType<T, Dim>::view_type ViewType;
        typedef typename detail::ViewType<T, Dim, Kokkos::MemoryTraits<Kokkos::Atomic>>::view_type
            AtomicViewType;
 
        // Constructors ///////////////////////////////////////////////////////
 
        LagrangeSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
                      const QuadratureType& quadrature, Layout_t& layout);
 
        LagrangeSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
                      const QuadratureType& quadrature);
 
        void initialize(UniformCartesian<T, Dim>& mesh, Layout_t& layout);
 
 
        void initializeElementIndices(Layout_t& layout);
 
 
        void updateLayout(Layout_t& layout);
 
 
        KOKKOS_FUNCTION size_t numGlobalDOFs() const override;
 
        KOKKOS_FUNCTION size_t getLocalDOFIndex(const size_t& elementIndex,
                                                const size_t& globalDOFIndex) const override;
 
        KOKKOS_FUNCTION size_t getGlobalDOFIndex(const size_t& elementIndex,
                                                 const size_t& localDOFIndex) const override;
 
        KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getLocalDOFIndices() const override;
 
        KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
            const size_t& element_index) const override;
 
 
        KOKKOS_FUNCTION T evaluateRefElementShapeFunction(const size_t& localDOF,
                                                          const point_t& localPoint) const;
 
        KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionGradient(
            const size_t& localDOF, const point_t& localPoint) const;
 
 
        KOKKOS_FUNCTION point_t
        getInverseTransposeTransformationJacobian(vertex_points_t pt) const {
            return this->ref_element_m.getInverseTransposeTransformationJacobian(pt);
        }
 
 
        template <typename F>
        FieldLHS evaluateAx(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_lower(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_upper(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_upperlower(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_inversediag(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_diag(FieldLHS& field, F& evalFunction);
 
        template <typename F>
        FieldLHS evaluateAx_lift(FieldLHS& field, F& evalFunction);
 
        void evaluateLoadVector(FieldRHS& field) const;
        void evaluateLumpedMass(FieldRHS& field) const;
 
 
        template <typename F>
        T computeErrorL2(const FieldLHS& u_h, const F& u_sol) const;
 
        T computeAvg(const FieldLHS& u_h) const;
 
        struct DeviceStruct {
            // members we need to copy for the following functions:
            // works since numElementDOFs in LagrangeSpace is static constexpr
            static constexpr unsigned numElementDOFs = LagrangeSpace::numElementDOFs;
            Vector<size_t, Dim> nr_m;
            ElementType ref_element_m;
 
            // these are the functions needed for interpolation to the space
            KOKKOS_FUNCTION indices_t getMeshVertexNDIndex(const size_t& vertex_index) const;
 
            KOKKOS_FUNCTION size_t getLocalDOFIndex(const indices_t& elementNDIndex,
                                                    const size_t& globalDOFIndex) const;
            KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
                const indices_t& elementNDIndex) const;
 
            KOKKOS_FUNCTION T evaluateRefElementShapeFunction(const size_t& localDOF,
                                                              const point_t& localPoint) const;
            KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionGradient(
                const size_t& localDOF, const point_t& localPoint) const;
        };
 
        DeviceStruct getDeviceMirror() const;
 
    private:
        KOKKOS_FUNCTION bool isDOFOnBoundary(const indices_t& ndindex) const {
            for (size_t d = 0; d < Dim; ++d) {
                if (ndindex[d] <= 0 || ndindex[d] >= this->nr_m[d] - 1) {
                    return true;
                }
            }
            return false;
        }
 
        Kokkos::View<size_t*> elementIndices;
 
        // One time allocated field of type FieldLHS to store results
        FieldLHS resultField;
    };
 
}  // namespace ippl
 
#include "FEM/LagrangeSpace.hpp"
 
#endif