NedelecSpace.h

// Class NedelecSpace
//    This is the NedelecSpace class. It is a class representing a Nedelec space
//    for finite element methods on a structured grid.
 
#ifndef IPPL_NEDELECSPACE_H
#define IPPL_NEDELECSPACE_H
 
#include <cmath>
 
#include "FEM/FEMQuadratureData.h"
#include "FEM/FEMVector.h"
#include "FEM/FiniteElementSpace.h"
 
constexpr unsigned getNedelecNumElementDOFs(unsigned Dim, [[maybe_unused]] unsigned Order) {
    // needs to be constexpr pow function to work at compile time. Kokkos::pow
    // doesn't work.
    return static_cast<unsigned>(static_cast<int>(Dim) * power(2, static_cast<int>(Dim - 1)));
}
 
namespace ippl {
 
    template <typename T, unsigned Dim, unsigned Order, typename ElementType,
              typename QuadratureType, typename FieldType>
    // requires IsQuadrature<QuadratureType>
    class NedelecSpace
        : public FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order), ElementType,
                                    QuadratureType, FEMVector<T>, FEMVector<T>> {
    public:
        // The number of degrees of freedom per element
        static constexpr unsigned numElementDOFs = getNedelecNumElementDOFs(Dim, Order);
 
        // The dimension of the mesh
        static constexpr unsigned dim =
            FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector<T>,
                               FEMVector<T>>::dim;
 
        // The order of the Nedelec space
        static constexpr unsigned order = Order;
 
        // The number of mesh vertices per element
        static constexpr unsigned numElementVertices =
            FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector<T>,
                               FEMVector<T>>::numElementVertices;
 
        // A vector with the position of the element in the mesh in each dimension
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FEMVector<T>, FEMVector<T>>::indices_t indices_t;
 
        // A point in the global coordinate system
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FEMVector<T>, FEMVector<T>>::point_t point_t;
 
        typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
                                            FEMVector<T>, FEMVector<T>>::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, 1>::view_type ViewType;
        typedef typename detail::ViewType<T, 1, Kokkos::MemoryTraits<Kokkos::Atomic>>::view_type
            AtomicViewType;
 
        // Constructors ///////////////////////////////////////////////////////
 
        NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
                     const QuadratureType& quadrature, const Layout_t& layout);
 
        NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
                     const QuadratureType& quadrature);
 
        void initialize(UniformCartesian<T, Dim>& mesh, const Layout_t& layout);
 
 
        void initializeElementIndices(const Layout_t& layout);
 
        NDIndex<Dim> getLocalNDIndex() const { return layout_m.getLocalNDIndex(); }
 
 
        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& elementIndex) const override;
 
        KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
            const indices_t& elementIndex) const;
 
        KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getFEMVectorDOFIndices(
            const size_t& elementIndex, NDIndex<Dim> ldom) const;
 
        KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getFEMVectorDOFIndices(
            indices_t elementIndex, NDIndex<Dim> ldom) const;
 
        KOKKOS_FUNCTION point_t getLocalDOFPosition(size_t localDOFIndex) const;
 
 
        KOKKOS_FUNCTION point_t evaluateRefElementShapeFunction(const size_t& localDOF,
                                                                const point_t& localPoint) const;
 
        KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionCurl(
            const size_t& localDOF, const point_t& localPoint) const;
 
 
        template <typename F>
        FEMVector<T> evaluateAx(FEMVector<T>& x, F& evalFunction) const;
 
        FEMVector<T> evaluateLoadVector(const FEMVector<point_t>& f) const;
 
        template <typename F>
        FEMVector<T> evaluateLoadVectorFunctor(const F& f) const;
 
 
        FEMVector<T> createFEMVector() const;
 
        Kokkos::View<point_t*> reconstructToPoints(const Kokkos::View<point_t*>& positions,
                                                   const FEMVector<T>& coef) const;
 
 
        template <typename F>
        T computeError(const FEMVector<T>& u_h, const F& u_sol) const;
 
        KOKKOS_FUNCTION bool isDOFOnBoundary(const size_t& dofIdx) const;
 
        KOKKOS_FUNCTION int getBoundarySide(const size_t& dofIdx) const;
 
    private:
        FEMVector<T> createFEMVector2d() const;
 
        FEMVector<T> createFEMVector3d() const;
 
        Kokkos::View<size_t*> elementIndices;
 
        Vector<point_t, 12> localDofPositions_m;
 
        Layout_t layout_m;
    };
 
}  // namespace ippl
 
#include "FEM/NedelecSpace.hpp"
 
#endif