greville_interpolation_points.hpp

// Copyright (C) The DDC development team, see COPYRIGHT.md file
//
// SPDX-License-Identifier: MIT
 
#pragma once
 
#include <cstddef>
#include <type_traits>
#include <vector>
 
#include <ddc/ddc.hpp>
 
#include "spline_boundary_conditions.hpp"
 
namespace ddc {
 
template <class BSplines, ddc::BoundCond BcXmin, ddc::BoundCond BcXmax>
class GrevilleInterpolationPoints
{
    using tag_type = typename BSplines::tag_type;
 
    template <class Sampling>
    struct IntermediateUniformSampling
        : UniformPointSampling<typename Sampling::continuous_dimension_type>
    {
    };
 
    template <class Sampling>
    struct IntermediateNonUniformSampling
        : NonUniformPointSampling<typename Sampling::continuous_dimension_type>
    {
    };
 
    template <class Sampling, typename U = BSplines, class = std::enable_if_t<U::is_uniform()>>
    static auto uniform_greville_points()
    {
        using SamplingImpl = typename Sampling::template Impl<Sampling, Kokkos::HostSpace>;
 
        double constexpr shift = (BSplines::degree() % 2 == 0) ? 0.5 : 0.0;
        double dx
                = (ddc::discrete_space<BSplines>().rmax() - ddc::discrete_space<BSplines>().rmin())
                  / ddc::discrete_space<BSplines>().ncells();
        return SamplingImpl(
                ddc::Coordinate<tag_type>(ddc::discrete_space<BSplines>().rmin() + shift * dx),
                ddc::Coordinate<tag_type>(dx));
    }
 
    template <class Sampling, typename U = BSplines, class = std::enable_if_t<!U::is_uniform()>>
    static auto non_uniform_greville_points()
    {
        using SamplingImpl = typename Sampling::template Impl<Sampling, Kokkos::HostSpace>;
 
        int n_greville_points = ddc::discrete_space<BSplines>().nbasis();
        if constexpr (U::is_periodic()) {
            n_greville_points += 1;
        }
 
        std::vector<double> greville_points(n_greville_points);
        ddc::DiscreteDomain<BSplines> bspline_domain
                = ddc::discrete_space<BSplines>().full_domain().take_first(
                        ddc::DiscreteVector<BSplines>(ddc::discrete_space<BSplines>().nbasis()));
 
        ddc::for_each(bspline_domain, [&](ddc::DiscreteElement<BSplines> ib) {
            // Define the Greville points from the bspline knots
            greville_points[ib.uid()] = 0.0;
            for (std::size_t i(0); i < BSplines::degree(); ++i) {
                greville_points[ib.uid()]
                        += ddc::discrete_space<BSplines>().get_support_knot_n(ib, i + 1);
            }
            greville_points[ib.uid()] /= BSplines::degree();
        });
 
        std::vector<double> temp_knots(BSplines::degree());
        // Use periodicity to ensure all points are in the domain
        if constexpr (U::is_periodic()) {
            int npoints(0);
            // Count the number of interpolation points that need shifting to preserve the ordering
            while (greville_points[npoints] < ddc::discrete_space<BSplines>().rmin()) {
                temp_knots[npoints]
                        = greville_points[npoints] + ddc::discrete_space<BSplines>().length();
                npoints++;
            }
            // Shift the points
            for (std::size_t i = 0; i < ddc::discrete_space<BSplines>().nbasis() - npoints; ++i) {
                greville_points[i] = greville_points[i + npoints];
            }
            for (int i = 0; i < npoints; ++i) {
                greville_points[ddc::discrete_space<BSplines>().nbasis() - npoints + i]
                        = temp_knots[i];
            }
 
            // Save a periodic point to initialise the domain size
            greville_points[n_greville_points - 1]
                    = greville_points[0] + ddc::discrete_space<BSplines>().length();
        }
 
        return SamplingImpl(greville_points);
    }
 
    static constexpr int N_BE_MIN = n_boundary_equations(BcXmin, BSplines::degree());
    static constexpr int N_BE_MAX = n_boundary_equations(BcXmax, BSplines::degree());
    template <class U>
    static constexpr bool is_uniform_mesh_v
            = U::is_uniform() && ((N_BE_MIN != 0 && N_BE_MAX != 0) || U::is_periodic());
 
public:
    template <
            class Sampling,
            typename U = BSplines,
            std::enable_if_t<
                    is_uniform_mesh_v<U>,
                    bool> = true> // U must be in condition for SFINAE
    static auto get_sampling()
    {
        return uniform_greville_points<Sampling>();
    }
 
    template <
            class Sampling,
            typename U = BSplines,
            std::enable_if_t<
                    !is_uniform_mesh_v<U>,
                    bool> = true> // U must be in condition for SFINAE
    static auto get_sampling()
    {
        using SamplingImpl = typename Sampling::template Impl<Sampling, Kokkos::HostSpace>;
        if constexpr (U::is_uniform()) {
            using IntermediateSampling = IntermediateUniformSampling<Sampling>;
            auto points_wo_bcs = uniform_greville_points<IntermediateSampling>();
            int const n_break_points = ddc::discrete_space<BSplines>().ncells() + 1;
            int const npoints = ddc::discrete_space<BSplines>().nbasis() - N_BE_MIN - N_BE_MAX;
            std::vector<double> points_with_bcs(npoints);
 
            // Construct Greville-like points at the edge
            if constexpr (BcXmin == ddc::BoundCond::GREVILLE) {
                for (std::size_t i(0); i < BSplines::degree() / 2 + 1; ++i) {
                    points_with_bcs[i]
                            = (BSplines::degree() - i) * ddc::discrete_space<BSplines>().rmin();
                    for (std::size_t j(0); j < i; ++j) {
                        points_with_bcs[i] += ddc::discrete_space<BSplines>().get_support_knot_n(
                                ddc::DiscreteElement<BSplines>(i),
                                BSplines::degree() - j);
                    }
                    points_with_bcs[i] /= BSplines::degree();
                }
            } else {
                points_with_bcs[0]
                        = points_wo_bcs.coordinate(ddc::DiscreteElement<IntermediateSampling>(0));
            }
 
            int const n_start
                    = (BcXmin == ddc::BoundCond::GREVILLE) ? BSplines::degree() / 2 + 1 : 1;
            int const domain_size = n_break_points - 2;
            ddc::DiscreteDomain<IntermediateSampling> const
                    domain(ddc::DiscreteElement<IntermediateSampling>(1),
                           ddc::DiscreteVector<IntermediateSampling>(domain_size));
 
            // Copy central points
            ddc::for_each(domain, [&](auto ip) {
                points_with_bcs[ip.uid() + n_start - 1] = points_wo_bcs.coordinate(ip);
            });
 
            // Construct Greville-like points at the edge
            if constexpr (BcXmax == ddc::BoundCond::GREVILLE) {
                for (std::size_t i(0); i < BSplines::degree() / 2 + 1; ++i) {
                    points_with_bcs[npoints - 1 - i]
                            = (BSplines::degree() - i) * ddc::discrete_space<BSplines>().rmax();
                    for (std::size_t j(0); j < i; ++j) {
                        points_with_bcs[npoints - 1 - i]
                                += ddc::discrete_space<BSplines>().get_support_knot_n(
                                        ddc::DiscreteElement<BSplines>(
                                                ddc::discrete_space<BSplines>().nbasis() - 1 - i),
                                        j + 1);
                    }
                    points_with_bcs[npoints - 1 - i] /= BSplines::degree();
                }
            } else {
                points_with_bcs[npoints - 1]
                        = points_wo_bcs.coordinate(ddc::DiscreteElement<IntermediateSampling>(
                                ddc::discrete_space<BSplines>().ncells() - 1
                                + BSplines::degree() % 2));
            }
            return SamplingImpl(points_with_bcs);
        } else {
            using IntermediateSampling = IntermediateNonUniformSampling<Sampling>;
            if constexpr (N_BE_MIN == 0 && N_BE_MAX == 0) {
                return non_uniform_greville_points<Sampling>();
            } else {
                auto points_wo_bcs = non_uniform_greville_points<IntermediateSampling>();
                // All points are Greville points. Extract unnecessary points near the boundary
                std::vector<double> points_with_bcs(points_wo_bcs.size() - N_BE_MIN - N_BE_MAX);
                int constexpr n_start = N_BE_MIN;
 
                using length = ddc::DiscreteVector<IntermediateSampling>;
 
                ddc::DiscreteDomain<IntermediateSampling> const
                        domain(ddc::DiscreteElement<IntermediateSampling>(n_start),
                               length(points_with_bcs.size()));
 
                points_with_bcs[0] = points_wo_bcs.coordinate(domain.front());
                ddc::for_each(domain.remove(length(1), length(1)), [&](auto ip) {
                    points_with_bcs[ip.uid() - n_start] = points_wo_bcs.coordinate(ip);
                });
                points_with_bcs[points_with_bcs.size() - 1]
                        = points_wo_bcs.coordinate(domain.back());
 
                return SamplingImpl(points_with_bcs);
            }
        }
    }
 
    using interpolation_mesh_type = std::conditional_t<
            is_uniform_mesh_v<BSplines>,
            ddc::UniformPointSampling<tag_type>,
            ddc::NonUniformPointSampling<tag_type>>;
 
    template <class Sampling>
    static ddc::DiscreteDomain<Sampling> get_domain()
    {
        int const npoints = ddc::discrete_space<BSplines>().nbasis() - N_BE_MIN - N_BE_MAX;
        return ddc::DiscreteDomain<Sampling>(
                ddc::DiscreteElement<Sampling>(0),
                ddc::DiscreteVector<Sampling>(npoints));
    }
};
} // namespace ddc