128 lines
4.9 KiB
128 lines
4.9 KiB
#pragma once
|
|
|
|
#include <Eigen/Core>
|
|
#include <random>
|
|
|
|
namespace meshless {
|
|
|
|
template <typename T>
|
|
struct Pi {
|
|
static constexpr T value = T(3.14159265358979323846264338327950L); ///< Value of @f$\pi@f$.
|
|
};
|
|
|
|
static const double PI = Pi<double>::value; ///< Mathematical constant pi in double precision.
|
|
static const double INF = 1.0 / 0.0; ///< Infinite floating point value.
|
|
static const double NaN = 0.0 / 0.0; ///< Not-a-number floating point value.
|
|
|
|
template <typename scalar_t, int dim>
|
|
struct SphereDiscretization {
|
|
/**
|
|
* Construct a randomized discretization.
|
|
* @param radius Radius of the sphere.
|
|
* @param num_samples Number of points on the equator, implies nodal spacing `dp = 2*pi*r/n`.
|
|
* @param generator A random number generator.
|
|
* @return A vector of discretization points.
|
|
*/
|
|
template <typename generator_t>
|
|
static std::vector<Eigen::Matrix<scalar_t, dim, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> construct(
|
|
scalar_t radius, int num_samples, generator_t& generator) {
|
|
scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
|
|
std::uniform_real_distribution<scalar_t> distribution(0, Pi<scalar_t>::value);
|
|
scalar_t offset = distribution(generator);
|
|
std::vector<Eigen::Matrix<scalar_t, dim, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> result;
|
|
for(int i = 0; i < num_samples / 2; ++i) {
|
|
scalar_t phi = i * dphi + offset;
|
|
if(phi > Pi<scalar_t>::value) phi -= Pi<scalar_t>::value;
|
|
int slice_n = static_cast<int>(std::ceil(num_samples * std::sin(phi)));
|
|
if(slice_n == 0) continue;
|
|
auto slice = SphereDiscretization<scalar_t, dim - 1>::construct(
|
|
radius * std::sin(phi), slice_n, generator);
|
|
Eigen::Matrix<scalar_t, dim, 1> v;
|
|
for(const auto& p : slice) {
|
|
v[0] = radius * std::cos(phi);
|
|
v.template tail<dim - 1>() = p;
|
|
result.push_back(v);
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
/// Construct the discretization.
|
|
static std::vector<Eigen::Matrix<scalar_t, dim, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> construct(
|
|
scalar_t radius, int num_samples) {
|
|
scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
|
|
std::vector<Eigen::Matrix<scalar_t, dim, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> result;
|
|
for(int i = 0; i < num_samples / 2; ++i) {
|
|
scalar_t phi = i * dphi;
|
|
if(phi > Pi<scalar_t>::value) phi -= Pi<scalar_t>::value;
|
|
int slice_n = static_cast<int>(std::ceil(num_samples * std::sin(phi)));
|
|
if(slice_n == 0) continue;
|
|
auto slice = SphereDiscretization<scalar_t, dim - 1>::construct(
|
|
radius * std::sin(phi), slice_n);
|
|
Eigen::Matrix<scalar_t, dim, 1> v;
|
|
for(const auto& p : slice) {
|
|
v[0] = radius * std::cos(phi);
|
|
v.template tail<dim - 1>() = p;
|
|
result.push_back(v);
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
};
|
|
|
|
/// Two-dimensional base case of the discretization.
|
|
template <typename scalar_t>
|
|
struct SphereDiscretization<scalar_t, 2> {
|
|
/// Construct a randomized discretization.
|
|
template <typename generator_t>
|
|
static std::vector<Eigen::Matrix<scalar_t, 2, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> construct(
|
|
scalar_t radius, int num_samples, generator_t& generator) {
|
|
scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
|
|
std::uniform_real_distribution<scalar_t> distribution(0, 2 * Pi<scalar_t>::value);
|
|
scalar_t offset = distribution(generator);
|
|
std::vector<Eigen::Matrix<scalar_t, 2, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> result;
|
|
for(int i = 0; i < num_samples; ++i) {
|
|
scalar_t phi = i * dphi + offset;
|
|
result.emplace_back(radius * std::cos(phi), radius * std::sin(phi));
|
|
}
|
|
return result;
|
|
}
|
|
/// Construct the discretization.
|
|
static std::vector<Eigen::Matrix<scalar_t, 2, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> construct(
|
|
scalar_t radius, int num_samples) {
|
|
scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
|
|
std::vector<Eigen::Matrix<scalar_t, 2, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> result;
|
|
for(int i = 0; i < num_samples; ++i) {
|
|
scalar_t phi = i * dphi;
|
|
result.emplace_back(radius * std::cos(phi), radius * std::sin(phi));
|
|
}
|
|
return result;
|
|
}
|
|
};
|
|
|
|
/// One-dimensional base case of the discretization.
|
|
template <typename scalar_t>
|
|
struct SphereDiscretization<scalar_t, 1> {
|
|
template <typename generator_t>
|
|
/// Construct a randomized discretization.
|
|
static std::vector<Eigen::Matrix<scalar_t, 1, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 1, 1>>> construct(
|
|
scalar_t radius, int, generator_t&) {
|
|
return { Eigen::Matrix<scalar_t, 1, 1>(-radius), Eigen::Matrix<scalar_t, 1, 1>(radius) };
|
|
}
|
|
/// Construct the discretization.
|
|
static std::vector<Eigen::Matrix<scalar_t, 1, 1>,
|
|
Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 1, 1>>> construct(
|
|
scalar_t radius, int) {
|
|
return { Eigen::Matrix<scalar_t, 1, 1>(-radius), Eigen::Matrix<scalar_t, 1, 1>(radius) };
|
|
}
|
|
};
|
|
|
|
};
|