NedelecSpace.h

Go to the documentation of this file.

// 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/FiniteElementSpace.h"
#include "FEM/FEMVector.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);


        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