Skip to content

Commit

Permalink
Align UniformBSplines API to NonUniformBSplines API (#489)
Browse files Browse the repository at this point in the history 
* Remove get_knot and get_support_knot_n

* Update get_first_support_knot and get_last_support_knot to match bsplines_non_uniform. Add knot domain as well as break point domain

* Improve variable name

* Update docstring

* Remove use of get_support_knot_n

* Correct domain

* Correct uniform knot mesh initialisation

* Use init_ghosted to create domains

* Apply suggestions from code review

---------

Co-authored-by: Emily Bourne <[email protected]>
Co-authored-by: Thomas Padioleau <[email protected]>
  • Loading branch information
@EmilyBourne @tpadioleau
3 people committed Jun 10, 2024
1 parent 4607762 commit b8ccaf4
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 67 deletions.
8 changes: 4 additions & 4 deletions include/ddc/kernels/splines/bsplines_non_uniform.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase
using discrete_vector_type = DiscreteVector<DDim>;

private:
ddc::DiscreteDomain<knot_mesh_type> m_domain;
ddc::DiscreteDomain<knot_mesh_type> m_knot_domain;
ddc::DiscreteDomain<knot_mesh_type> m_break_point_domain;

public:
Expand Down Expand Up @@ -161,7 +161,7 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase
*/
template <class OriginMemorySpace>
explicit Impl(Impl<DDim, OriginMemorySpace> const& impl)
: m_domain(impl.m_domain)
: m_knot_domain(impl.m_knot_domain)
, m_break_point_domain(impl.m_break_point_domain)
{
}
Expand Down Expand Up @@ -352,7 +352,7 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase
*/
KOKKOS_INLINE_FUNCTION std::size_t npoints() const noexcept
{
return m_domain.size() - 2 * degree();
return m_knot_domain.size() - 2 * degree();
}

/** @brief Returns the number of basis functions.
Expand Down Expand Up @@ -414,7 +414,7 @@ template <class RandomIt>
NonUniformBSplines<Tag, D>::Impl<DDim, MemorySpace>::Impl(
RandomIt const break_begin,
RandomIt const break_end)
: m_domain(
: m_knot_domain(
ddc::DiscreteElement<knot_mesh_type>(0),
ddc::DiscreteVector<knot_mesh_type>(
(break_end - break_begin)
Expand Down
79 changes: 28 additions & 51 deletions include/ddc/kernels/splines/bsplines_uniform.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,8 +104,9 @@ class UniformBSplines : detail::UniformBSplinesBase
using discrete_vector_type = DiscreteVector<DDim>;

private:
// In the periodic case, it contains twice the periodic point!!!
ddc::DiscreteDomain<knot_mesh_type> m_domain;
// In the periodic case, they contain the periodic point twice!!!
ddc::DiscreteDomain<knot_mesh_type> m_knot_domain;
ddc::DiscreteDomain<knot_mesh_type> m_break_point_domain;

public:
Impl() = default;
Expand All @@ -117,24 +118,28 @@ class UniformBSplines : detail::UniformBSplinesBase
* @param ncells the number of cells in the range [rmin, rmax]
*/
explicit Impl(ddc::Coordinate<Tag> rmin, ddc::Coordinate<Tag> rmax, std::size_t ncells)
: m_domain(
ddc::DiscreteElement<knot_mesh_type>(0),
ddc::DiscreteVector<knot_mesh_type>(
ncells + 1)) // Create a mesh including the eventual periodic point
{
assert(ncells > 0);
ddc::init_discrete_space<knot_mesh_type>(knot_mesh_type::template init<knot_mesh_type>(
rmin,
rmax,
ddc::DiscreteVector<knot_mesh_type>(ncells + 1)));
ddc::DiscreteDomain<knot_mesh_type> pre_ghost;
ddc::DiscreteDomain<knot_mesh_type> post_ghost;
std::tie(m_break_point_domain, m_knot_domain, pre_ghost, post_ghost)
= ddc::init_discrete_space<knot_mesh_type>(
knot_mesh_type::template init_ghosted<knot_mesh_type>(
rmin,
rmax,
ddc::DiscreteVector<knot_mesh_type>(ncells + 1),
ddc::DiscreteVector<knot_mesh_type>(degree()),
ddc::DiscreteVector<knot_mesh_type>(degree())));
}

/** @brief Copy-constructs from another Impl with a different Kokkos memory space
*
* @param impl A reference to the other Impl
*/
template <class OriginMemorySpace>
explicit Impl(Impl<DDim, OriginMemorySpace> const& impl) : m_domain(impl.m_domain)
explicit Impl(Impl<DDim, OriginMemorySpace> const& impl)
: m_knot_domain(impl.m_knot_domain)
, m_break_point_domain(impl.m_break_point_domain)
{
}

Expand Down Expand Up @@ -232,33 +237,19 @@ class UniformBSplines : detail::UniformBSplinesBase
integrals(ddc::ChunkSpan<double, discrete_domain_type, Layout, MemorySpace2>
int_vals) const;

/** @brief Returns the coordinate of the knot corresponding to the given index.
*
* Returns the coordinate of the knot corresponding to the given index. The domain
* over which the B-splines are defined is comprised of ncells+1 break points however there are a total of
* ncells+1+2*degree knots. The additional knots which control the shape of the B-splines near the
* boundary are added equidistantly before and after the break points. The knot index is therefore in the interval [-degree, ncells+degree]
*
* @param[in] knot_idx Integer identifying index of the knot.
* @return Coordinate of the knot.
*/
KOKKOS_INLINE_FUNCTION ddc::Coordinate<Tag> get_knot(int knot_idx) const noexcept
{
return ddc::Coordinate<Tag>(rmin() + knot_idx * ddc::step<knot_mesh_type>());
}

/** @brief Returns the coordinate of the first support knot associated to a DiscreteElement identifying a B-spline.
*
* Each B-spline has a support defined over (degree+2) knots. For a B-spline identified by the
* provided DiscreteElement, this function returns the first knot in the support of the B-spline.
* In other words it returns the lower bound of the support.
*
* @param[in] ix DiscreteElement identifying the B-spline.
* @return Coordinate of the knot.
* @return DiscreteElement of the lower bound of the support of the B-spline.
*/
KOKKOS_INLINE_FUNCTION double get_first_support_knot(discrete_element_type const& ix) const
KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<knot_mesh_type> get_first_support_knot(
discrete_element_type const& ix) const
{
return get_knot(ix.uid() - degree());
return ddc::DiscreteElement<knot_mesh_type>((ix - discrete_element_type(0)).value());
}

/** @brief Returns the coordinate of the last support knot associated to a DiscreteElement identifying a B-spline.
Expand All @@ -268,26 +259,12 @@ class UniformBSplines : detail::UniformBSplinesBase
* In other words it returns the upper bound of the support.
*
* @param[in] ix DiscreteElement identifying the B-spline.
* @return Coordinate of the knot.
*/
KOKKOS_INLINE_FUNCTION double get_last_support_knot(discrete_element_type const& ix) const
{
return get_knot(ix.uid() + 1);
}

/** @brief Returns the coordinate of the (n+1)-th knot in the support of the identified B-spline.
*
* Each B-spline has a support defined over (degree+2) knots. For a B-spline identified by the
* provided DiscreteElement, this function returns the (n+1)-th knot in the support of the B-spline.
*
* @param[in] ix DiscreteElement identifying the B-spline.
* @param[in] n Integer indexing a knot in the support of the B-spline.
* @return Coordinate of the knot.
* @return DiscreteElement of the upper bound of the support of the B-spline.
*/
KOKKOS_INLINE_FUNCTION double get_support_knot_n(discrete_element_type const& ix, int n)
const
KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<knot_mesh_type> get_last_support_knot(
discrete_element_type const& ix) const
{
return get_knot(ix.uid() + n - degree());
return get_first_support_knot(ix) + ddc::DiscreteVector<knot_mesh_type>(degree() + 1);
}

/** @brief Returns the coordinate of the lower bound of the domain on which the B-splines are defined.
Expand All @@ -296,7 +273,7 @@ class UniformBSplines : detail::UniformBSplinesBase
*/
KOKKOS_INLINE_FUNCTION ddc::Coordinate<Tag> rmin() const noexcept
{
return ddc::coordinate(m_domain.front());
return ddc::coordinate(m_break_point_domain.front());
}

/** @brief Returns the coordinate of the upper bound of the domain on which the B-splines are defined.
Expand All @@ -305,7 +282,7 @@ class UniformBSplines : detail::UniformBSplinesBase
*/
KOKKOS_INLINE_FUNCTION ddc::Coordinate<Tag> rmax() const noexcept
{
return ddc::coordinate(m_domain.back());
return ddc::coordinate(m_break_point_domain.back());
}

/** @brief Returns the length of the domain.
Expand Down Expand Up @@ -345,7 +322,7 @@ class UniformBSplines : detail::UniformBSplinesBase
*/
KOKKOS_INLINE_FUNCTION ddc::DiscreteDomain<knot_mesh_type> break_point_domain() const
{
return m_domain;
return m_break_point_domain;
}

/** @brief Returns the number of basis functions.
Expand All @@ -368,7 +345,7 @@ class UniformBSplines : detail::UniformBSplinesBase
*/
KOKKOS_INLINE_FUNCTION std::size_t ncells() const noexcept
{
return m_domain.size() - 1;
return m_break_point_domain.size() - 1;
}

private:
Expand Down
34 changes: 22 additions & 12 deletions include/ddc/kernels/splines/greville_interpolation_points.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,11 +164,17 @@ class GrevilleInterpolationPoints
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);
}
ddc::DiscreteElement<BSplines> spline_idx(i);
ddc::DiscreteVector<UniformBsplinesKnots<BSplines>> n_knots_in_domain(i);
ddc::DiscreteDomain<UniformBsplinesKnots<BSplines>> sub_domain(
ddc::discrete_space<BSplines>().get_last_support_knot(spline_idx)
- n_knots_in_domain,
n_knots_in_domain);
ddc::for_each(
sub_domain,
[&](DiscreteElement<UniformBsplinesKnots<BSplines>> ik) {
points_with_bcs[i] += ddc::coordinate(ik);
});
points_with_bcs[i] /= BSplines::degree();
}
} else {
Expand All @@ -193,13 +199,17 @@ class GrevilleInterpolationPoints
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);
}
ddc::DiscreteElement<BSplines> spline_idx(
ddc::discrete_space<BSplines>().nbasis() - 1 - i);
ddc::DiscreteVector<UniformBsplinesKnots<BSplines>> n_knots_in_domain(i);
ddc::DiscreteDomain<UniformBsplinesKnots<BSplines>> sub_domain(
ddc::discrete_space<BSplines>().get_first_support_knot(spline_idx) + 1,
n_knots_in_domain);
ddc::for_each(
sub_domain,
[&](DiscreteElement<UniformBsplinesKnots<BSplines>> ik) {
points_with_bcs[npoints - 1 - i] += ddc::coordinate(ik);
});
points_with_bcs[npoints - 1 - i] /= BSplines::degree();
}
} else {
Expand Down

0 comments on commit b8ccaf4

Please sign in to comment.