You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
307 lines
13 KiB
307 lines
13 KiB
#ifndef MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
|
|
#define MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
|
|
|
|
/**
|
|
* @file
|
|
* Implementation of NURBSShape class.
|
|
*/
|
|
|
|
#include "NURBSShape_fwd.hpp"
|
|
#include "GeneralSurfaceFill.hpp"
|
|
|
|
#include <medusa/bits/domains/DomainDiscretization.hpp>
|
|
#include <medusa/bits/utils/assert.hpp>
|
|
#include <medusa/bits/spatial_search/KDTreeMutable.hpp>
|
|
#include <medusa/bits/domains/NURBSPatch.hpp>
|
|
#include <medusa/bits/domains/BoxShape.hpp>
|
|
#include <medusa/bits/domains/compute_normal.hpp>
|
|
|
|
namespace mm
|
|
{
|
|
|
|
template <typename vec_t, typename param_vec_t>
|
|
NURBSShape<vec_t, param_vec_t>::NURBSShape(const Range<NURBSPatch<vec_t, param_vec_t>> &patches_in)
|
|
: patches{patches_in}, seed_(get_seed())
|
|
{
|
|
if (param_dim == 1)
|
|
{
|
|
n_samples = 2;
|
|
}
|
|
for (NURBSPatch<vec_t, param_vec_t> &patch : patches)
|
|
{
|
|
patch.computeDerivativeStructure();
|
|
}
|
|
}
|
|
|
|
template <typename vec_t, typename param_vec_t>
|
|
DomainDiscretization<vec_t> NURBSShape<vec_t, param_vec_t>::discretizeBoundaryWithDensity(
|
|
const std::function<scalar_t(vec_t)> &h, int type) const
|
|
{
|
|
return nurbs_shape_internal::NURBSShapeHelper<vec_t, param_vec_t>::
|
|
discretizeBoundaryWithDensity(*this, h, type);
|
|
}
|
|
|
|
template <typename vec_t, typename param_vec_t>
|
|
NURBSShape<vec_t, param_vec_t> &NURBSShape<vec_t, param_vec_t>::proximityTolerance(
|
|
scalar_t zeta)
|
|
{
|
|
assert_msg((0 < zeta && zeta < 1), "Zeta must be between 0 and 1, got %f.", zeta);
|
|
this->zeta = zeta;
|
|
return *this;
|
|
}
|
|
|
|
template <typename vec_t, typename param_vec_t>
|
|
NURBSShape<vec_t, param_vec_t> &NURBSShape<vec_t, param_vec_t>::boundaryProximity(
|
|
scalar_t epsilon)
|
|
{
|
|
assert_msg((0 < epsilon && epsilon < 1), "Epsilon must be between 0 and 1, got %f.", epsilon);
|
|
this->epsilon = epsilon;
|
|
return *this;
|
|
}
|
|
|
|
/// @cond
|
|
template <typename vec_t>
|
|
struct nurbs_shape_internal::NURBSShapeHelper<vec_t, Vec<typename vec_t::scalar_t, 2>>
|
|
{
|
|
typedef Vec<typename vec_t::scalar_t, 2> param_vec_t;
|
|
typedef typename vec_t::scalar_t scalar_t;
|
|
|
|
enum
|
|
{
|
|
dim = vec_t::dim,
|
|
proj_dim = dim + 1,
|
|
param_dim = param_vec_t::dim
|
|
};
|
|
typedef Vec<scalar_t, proj_dim> proj_vec_t;
|
|
|
|
static DomainDiscretization<vec_t> discretizeBoundaryWithDensity(
|
|
const NURBSShape<vec_t, param_vec_t> &shape,
|
|
const std::function<scalar_t(vec_t)> &h, int type)
|
|
{
|
|
if (type == 0)
|
|
type = -1;
|
|
std::mt19937 gen(shape.seed_);
|
|
|
|
GeneralSurfaceFill<vec_t, param_vec_t> gsf;
|
|
gsf.seed(shape.seed_).maxPoints(shape.max_points).proximityTolerance(shape.zeta).numSamples(shape.n_samples);
|
|
|
|
// Global domain and KDTree.
|
|
DomainDiscretization<vec_t> domain(shape);
|
|
KDTreeMutable<vec_t> tree_global;
|
|
|
|
int cnt = 0;
|
|
// Discretize every patch.
|
|
for (const NURBSPatch<vec_t, param_vec_t> &patch : shape.patches)
|
|
{
|
|
std::cout << "正在边界离散化第：" << ++cnt << "张曲面\n";
|
|
BoxShape<param_vec_t> bs = patch.getDomain();
|
|
DomainDiscretization<param_vec_t> local_param_d(bs);
|
|
|
|
// Discretize boundaries first.
|
|
Range<NURBSPatch<vec_t, Vec<scalar_t, 1>>> boundaries = patch.getBoundaries();
|
|
int i = -1;
|
|
for (NURBSPatch<vec_t, Vec<scalar_t, 1>> &boundary : boundaries)
|
|
{
|
|
i++;
|
|
boundary.computeDerivativeStructure();
|
|
|
|
KDTreeMutable<vec_t> tree_local; // Local KD-Tree should be per-boundary.
|
|
// Generate seed node.
|
|
BoxShape<Vec<scalar_t, 1>> boundary_bs = boundary.getDomain();
|
|
Vec<scalar_t, 1> t_min = boundary_bs.beg();
|
|
Vec<scalar_t, 1> t_max = boundary_bs.end();
|
|
|
|
std::uniform_real_distribution<> dis(t_min(0), t_max(0));
|
|
Vec<scalar_t, 1> t_seed = Vec<scalar_t, 1>(dis(gen));
|
|
|
|
// Insert into local structures.
|
|
param_vec_t param_seed = patch.getPatchParameterFromBoundaryParameter(t_seed,
|
|
i, shape.epsilon);
|
|
scalar_t w_seed;
|
|
// std::cout << "t-Range & t: " << t_min(0) << "," << t_max(0) << ", " << t_seed(0) << '\n';
|
|
vec_t pt_seed = boundary.evaluate(t_seed, &w_seed);
|
|
local_param_d.addInternalNode(param_seed, 1);
|
|
tree_local.insert(pt_seed);
|
|
|
|
// Insert into global structures.
|
|
scalar_t check_radius_seed = h(pt_seed);
|
|
scalar_t d_sq = tree_global.size() == 0 ? 10 * check_radius_seed * check_radius_seed
|
|
: tree_global.query(pt_seed).second[0];
|
|
if (d_sq >= (shape.zeta * check_radius_seed) * (shape.zeta * check_radius_seed))
|
|
{
|
|
vec_t normal_seed = surface_fill_internal::
|
|
compute_normal(patch.jacobian(param_seed, pt_seed, w_seed));
|
|
domain.addBoundaryNode(pt_seed, type, normal_seed);
|
|
tree_global.insert(pt_seed);
|
|
}
|
|
|
|
// General surface fill the boundary.
|
|
int cur_node = local_param_d.size() - 1;
|
|
int end_node = local_param_d.size();
|
|
while (cur_node < end_node && end_node < shape.max_points)
|
|
{
|
|
param_vec_t param = local_param_d.pos(cur_node);
|
|
Vec<scalar_t, 1> t = patch.getBoundaryParameterFromPatchParameter(param, i);
|
|
vec_t pt, der;
|
|
std::tie(pt, der) = boundary.evaluatePointAndJacobian(t);
|
|
|
|
std::array<Vec<scalar_t, 1>, 2> candidates{Vec<scalar_t, 1>(1),
|
|
Vec<scalar_t, 1>(-1)};
|
|
|
|
// Filter candidates.
|
|
scalar_t alpha = h(pt) / der.norm();
|
|
for (const auto &u_cand : candidates)
|
|
{
|
|
Vec<scalar_t, 1> t_new = t + alpha * u_cand; // Generate candidate.
|
|
param_vec_t param_new = patch.getPatchParameterFromBoundaryParameter(
|
|
t_new, i, shape.epsilon);
|
|
if (!local_param_d.contains(param_new))
|
|
continue;
|
|
|
|
scalar_t w_new;
|
|
vec_t pt_new = boundary.evaluate(t_new, &w_new);
|
|
|
|
// Check radius must be new radius.
|
|
scalar_t check_radius = (pt - pt_new).norm();
|
|
|
|
// Check for local insertion.
|
|
scalar_t d_sq_loc = tree_local.query(pt_new).second[0];
|
|
if (d_sq_loc >= (shape.zeta * check_radius) *
|
|
(shape.zeta * check_radius))
|
|
{
|
|
local_param_d.addInternalNode(param_new, 1);
|
|
tree_local.insert(pt_new);
|
|
end_node++;
|
|
|
|
// Check for global insertion.
|
|
scalar_t d_sq_global = tree_global.query(pt_new).second[0];
|
|
if (d_sq_global >= (shape.zeta * check_radius) *
|
|
(shape.zeta * check_radius))
|
|
{
|
|
vec_t normal_new = surface_fill_internal::
|
|
compute_normal(patch.jacobian(param_new, pt_new, w_new));
|
|
domain.addBoundaryNode(pt_new, type, normal_new);
|
|
tree_global.insert(pt_new);
|
|
}
|
|
}
|
|
}
|
|
cur_node++;
|
|
}
|
|
}
|
|
|
|
// Random another parameter seed around the middle of the domain, so that fill is
|
|
// successful even if patches overlap.
|
|
std::uniform_real_distribution<> dis(0.0, 1.0);
|
|
param_vec_t param_another = bs.beg() + 0.25 * (bs.end() - bs.beg()) +
|
|
dis(gen) * 0.25 * (bs.end() - bs.beg());
|
|
local_param_d.addInternalNode(param_another, 1);
|
|
|
|
scalar_t w_another;
|
|
// std::cout << "start anao\n";
|
|
vec_t pt_another = patch.evaluate(param_another, &w_another);
|
|
scalar_t check_radius_another = h(pt_another);
|
|
scalar_t d_sq_another = tree_global.query(pt_another).second[0];
|
|
if (d_sq_another >= (shape.zeta * check_radius_another) *
|
|
(shape.zeta * check_radius_another))
|
|
{
|
|
vec_t normal_another = surface_fill_internal::
|
|
compute_normal(patch.jacobian(param_another, pt_another, w_another));
|
|
domain.addBoundaryNode(pt_another, type, normal_another);
|
|
tree_global.insert(pt_another);
|
|
}
|
|
// std::cout << "end anao\n";
|
|
|
|
// Surface fill.
|
|
auto param = [&patch, &bs](const param_vec_t &t)
|
|
{
|
|
// std::cout << bs.beg() << '\n';
|
|
// std::cout << bs.end() << '\n';
|
|
// std::cout << "in lamba: " << t << '\n';
|
|
return patch.evaluatePointAndJacobian(t);
|
|
};
|
|
// std::cout << "start fill\n";
|
|
gsf.fillParametrization(domain, local_param_d, param, h, tree_global, type);
|
|
// std::cout << "end fill\n";
|
|
}
|
|
|
|
return domain;
|
|
}
|
|
};
|
|
|
|
template <typename vec_t>
|
|
struct nurbs_shape_internal::NURBSShapeHelper<vec_t, Vec<typename vec_t::scalar_t, 1>>
|
|
{
|
|
typedef Vec<typename vec_t::scalar_t, 1> param_vec_t;
|
|
typedef typename vec_t::scalar_t scalar_t;
|
|
enum
|
|
{
|
|
dim = vec_t::dim,
|
|
proj_dim = dim + 1,
|
|
param_dim = param_vec_t::dim
|
|
};
|
|
typedef Vec<scalar_t, proj_dim> proj_vec_t;
|
|
|
|
static DomainDiscretization<vec_t> discretizeBoundaryWithDensity(
|
|
const NURBSShape<vec_t, param_vec_t> &shape,
|
|
const std::function<scalar_t(vec_t)> &h, int type)
|
|
{
|
|
if (type == 0)
|
|
type = -1;
|
|
std::mt19937 gen(shape.seed_);
|
|
|
|
GeneralSurfaceFill<vec_t, param_vec_t> gsf;
|
|
gsf.seed(shape.seed_).maxPoints(shape.max_points).proximityTolerance(shape.zeta).numSamples(shape.n_samples);
|
|
|
|
// Global domain and KDTree.
|
|
DomainDiscretization<vec_t> domain(shape);
|
|
KDTreeMutable<vec_t> tree_global;
|
|
|
|
// Discretize every patch.
|
|
for (const NURBSPatch<vec_t, param_vec_t> &patch : shape.patches)
|
|
{
|
|
BoxShape<param_vec_t> bs = patch.getDomain();
|
|
DomainDiscretization<param_vec_t> local_param_d(bs);
|
|
|
|
// Insert boundary edges and random parameter as seeds.
|
|
std::array<param_vec_t, 2> param_seeds{
|
|
patch.getPatchParameterFromBoundaryParameter({}, 0, shape.epsilon),
|
|
patch.getPatchParameterFromBoundaryParameter({}, 1, shape.epsilon)};
|
|
// Another random seed could be generated here to make this algorithm also work on
|
|
// overlapping curves. This would, however, create two "gaps" of size between h
|
|
// and 2h where the fronts meet, instead of just one.
|
|
|
|
for (const Vec1d ¶m_seed : param_seeds)
|
|
{
|
|
scalar_t w_seed;
|
|
vec_t pt_seed = patch.evaluate(param_seed, &w_seed);
|
|
local_param_d.addInternalNode(param_seed, 1);
|
|
|
|
// Insert into global structure.
|
|
scalar_t check_radius_seed = h(pt_seed);
|
|
scalar_t d_sq = tree_global.size() == 0 ? 10 * check_radius_seed * check_radius_seed
|
|
: tree_global.query(pt_seed).second[0];
|
|
if (d_sq >= (shape.zeta * check_radius_seed) * (shape.zeta * check_radius_seed))
|
|
{
|
|
vec_t normal_seed = surface_fill_internal::
|
|
compute_normal(patch.jacobian(param_seed, pt_seed, w_seed));
|
|
domain.addBoundaryNode(pt_seed, type, normal_seed);
|
|
tree_global.insert(pt_seed);
|
|
}
|
|
}
|
|
|
|
// Surface fill.
|
|
auto param = [&patch](const param_vec_t &t)
|
|
{
|
|
return patch.evaluatePointAndJacobian(t);
|
|
};
|
|
gsf.fillParametrization(domain, local_param_d, param, h, tree_global, type);
|
|
}
|
|
|
|
return domain;
|
|
}
|
|
};
|
|
/// @endcond
|
|
|
|
} // namespace mm
|
|
|
|
#endif // MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
|
|
|