FEMVector.h

// Class FEMVector
//    This class represents a one dimensional vector which can be used in the 
//    context of FEM to represent a field defined on the DOFs of a mesh.
 
 
#ifndef IPPL_FEMVECTOR_H
#define IPPL_FEMVECTOR_H
 
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "Types/ViewTypes.h"
#include "Field/HaloCells.h"
 
namespace ippl {
 
    template <typename T>
    class FEMVector : public detail::Expression<
                        FEMVector<T>,
                        sizeof(typename detail::ViewType<T, 1>::view_type)>{
    public:
        static constexpr unsigned dim = 1;
        
        using value_type = T;
 
        FEMVector(size_t n, std::vector<size_t> neighbors,
            std::vector< Kokkos::View<size_t*> > sendIdxs,
            std::vector< Kokkos::View<size_t*> > recvIdxs);
        
        
        KOKKOS_FUNCTION FEMVector(const FEMVector<T>& other);
 
        
        FEMVector(size_t n);
        
 
        void fillHalo();
 
        void accumulateHalo();
 
        void setHalo(T setValue);
 
 
        FEMVector<T>& operator= (T value);
 
 
        template <typename E, size_t N>
        FEMVector<T>& operator= (const detail::Expression<E, N>& expr);
 
        
        FEMVector<T>& operator= (const FEMVector<T>& v);
 
 
 
        KOKKOS_INLINE_FUNCTION T operator[] (size_t i) const;
        
 
        KOKKOS_INLINE_FUNCTION T operator() (size_t i) const;
 
 
        const Kokkos::View<T*>& getView() const;
 
 
        size_t size() const;
 
 
        FEMVector<T> deepCopy() const;
        
 
        template <typename K>
        FEMVector<K> skeletonCopy() const;
        
 
        void pack(const Kokkos::View<size_t*>& idxStore);
        
        
        template <typename Op>
        void unpack(const Kokkos::View<size_t*>& idxStore);
 
 
        struct Assign{
            KOKKOS_INLINE_FUNCTION void operator()(T& lhs, const T& rhs) const {
                lhs = rhs;
            }   
        };
 
 
        struct AssignAdd{
            KOKKOS_INLINE_FUNCTION void operator()(T& lhs, const T& rhs) const {
                lhs += rhs;
            }   
        };
 
 
    private:
        struct BoundaryInfo {
 
            BoundaryInfo (std::vector<size_t> neighbors,
                std::vector< Kokkos::View<size_t*> > sendIdxs,
                std::vector< Kokkos::View<size_t*> > recvIdxs_m);
 
            
            std::vector<size_t> neighbors_m;
 
            std::vector< Kokkos::View<size_t*> > sendIdxs_m;
 
            std::vector< Kokkos::View<size_t*> > recvIdxs_m;
 
 
            detail::FieldBufferData<T> commBuffer_m;
        };
        
        Kokkos::View<T*> data_m;
 
 
        std::shared_ptr<BoundaryInfo> boundaryInfo_m;
 
    };
 
 
    template <typename T>
    T innerProduct(const FEMVector<T>& a, const FEMVector<T>& b) {
        T localSum = 0;
        auto aView = a.getView();
        auto bView = b.getView();
 
        
        size_t n = aView.extent(0);
        Kokkos::parallel_reduce("FEMVector innerProduct", n,
            KOKKOS_LAMBDA(const size_t i, T& val){
                val += aView(i)*bView(i);
            },
            Kokkos::Sum<T>(localSum)
        );
 
        T globalSum = 0;
        ippl::Comm->allreduce(localSum, globalSum, 1, std::plus<T>());
        return globalSum;
    }
 
    template <typename T>
    T norm(const FEMVector<T>& v, int p = 2) {
 
        T local = 0;
        auto view = v.getView();
        size_t n = view.extent(0);
        switch (p) {
            case 0: {
                Kokkos::parallel_reduce("FEMVector l0 norm", n,
                    KOKKOS_LAMBDA(const size_t i, T& val) {
                        val = Kokkos::max(val, Kokkos::abs(view(i)));
                    },
                    Kokkos::Max<T>(local)
                );
                T globalMax = 0;
                ippl::Comm->allreduce(local, globalMax, 1, std::greater<T>());
                return globalMax;
            }
            default: {
                Kokkos::parallel_reduce("FEMVector lp norm", n,
                    KOKKOS_LAMBDA(const size_t i, T& val) {
                        val += Kokkos::pow(Kokkos::abs(view(i)), p);
                    },
                    Kokkos::Sum<T>(local)
                );
                T globalSum = 0;
                ippl::Comm->allreduce(local, globalSum, 1, std::plus<T>());
                return std::pow(globalSum, 1.0 / p);
            }
        }
    }
}   // namespace ippl
 
 
#include "FEMVector.hpp"
 
#endif  // IPPL_FEMVECTOR_H