blender
/
blender
136
407
Fork 589
Code Issues 6.3k Pull Requests 768 Projects 19 Wiki Activity

Curves: initial surface collision for curves sculpt mode #104469

Merged
Jacques Lucke merged 29 commits from JacquesLucke/blender:temp-curves-surface-collision into main 2023-02-11 13:46:39 +01:00
Conversation 3 Commits 29 Files Changed 13 -1126
7 changed files with 0 additions and 1126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit e8991adbfb - Show all commits

2
source/blender/editors/sculpt_paint/curves_sculpt_pinch.cc
View File

@ -23,8 +23,6 @@
#include "DNA_screen_types.h"
#include "DNA_space_types.h"
#include "GEO_constraint_solver.hh"
#include "ED_screen.h"
#include "ED_view3d.h"

17
source/blender/geometry/CMakeLists.txt
View File

@ -16,7 +16,6 @@ set(INC
set(SRC
intern/add_curves_on_mesh.cc
intern/constraint_solver.cc
intern/curve_constraint_solver.cc
intern/fillet_curves.cc
intern/mesh_merge_by_distance.cc
@ -34,7 +33,6 @@ set(SRC
intern/uv_parametrizer.cc
GEO_add_curves_on_mesh.hh
GEO_constraint_solver.hh
GEO_curve_constraint_solver.hh
GEO_fillet_curves.hh
GEO_mesh_merge_by_distance.hh
@ -84,18 +82,3 @@ if(WITH_TBB)
endif()
blender_add_lib(bf_geometry "${SRC}" "${INC}" "${INC_SYS}" "${LIB}")
if(WITH_GTESTS)
set(TEST_SRC
tests/GEO_constraint_solver_test.cc
)
set(TEST_INC
)
set(TEST_LIB
bf_geometry
)
include(GTestTesting)
blender_add_test_executable(geometry "${TEST_SRC}" "${INC};${TEST_INC}" "${INC_SYS}" "${LIB};${TEST_LIB}")
add_subdirectory(tests/performance)
endif()

188
source/blender/geometry/GEO_constraint_solver.hh
View File

@ -1,188 +0,0 @@
/* SPDX-License-Identifier: GPL-2.0-or-later */
#pragma once
#include "BKE_bvhutils.h"
#include "BKE_curves.hh"
/** \file
* \ingroup geo
* \brief Constraint solver for curve deformations.
*/
namespace blender::geometry {
class ConstraintSolver {
public:
enum SolverType {
/* A fast single-iteration solver that solves constraints sequentially
* from the root down (Follow-the-Leader, FTL).
* The result is not physically correct, but very fast to compute.
* This solver type is suitable for editing, but not so much for accurate
* physical simulations, since the root side of a curve is infinitely stiff.
*
* For details see "Fast Simulation of Inextensible Hair and Fur"
* by Matthias Mueller and Tae Yong Kim. */
Sequential,
};
struct Params {
SolverType solver_type = SolverType::Sequential;
/* Keep the distance between points constant. */
bool use_length_constraints = true;
/* Root point is fixed to the surface. */
bool use_root_constraints = true;
/* Points do not penetrate the surface. */
bool use_collision_constraints = true;
/* Compliance (inverse stiffness)
* Alpha is used in physical simulation to control the softness of a constraint:
* For alpha == 0 the constraint is stiff and the maximum correction factor is applied.
* For values > 0 the constraint becomes squishy, and some violation is
* permitted, and the constraint gets corrected over multiple time steps.
*/
float alpha = 0.0f;
/* Number of substeps to perform.
* More substeps can be faster overall because of reduced search radius for collisions. */
int substep_count = 20;
/* Maximum overall distance a particle can move during a step.
* Divide by substep count to get max substep travel.
* This determines a the search radius for collisions.
* A larger travel distance means the point can move faster,
* but it can take longer to find collisions. */
float max_travel_distance = 0.1f;
/* Maximum number of simultaneous contacts to record per point. */
int max_contacts_per_point = 4;
/* Branching factor for the surface mesh BVH. */
int bvh_branching_factor = 2;
/* Number of iterations to satisfy constraints. */
int max_solver_iterations = 5;
/* Acceptable error threshold for convergence, in length units. */
float error_threshold = 1.0e-4f;
};
struct Result {
enum Status {
Ok,
ErrorNoConvergence,
};
Status status = Status::Ok;
/* True if any point's original travel was larger than the allowed maximum.
* If clamping is enabled the point's travel will be shorter than the input. */
bool max_travel_exceeded = false;
/* Total number of constraints solved. */
int constraint_count = 0;
/* Residual error values to judge convergence.
* Eror computation must be explicitly enabled during solving. */
struct {
/* Root-mean-square of the constraint residuals, indicating numerical quality
* of the solution. For a positional solver this value is in length units
* and describes how far constraints are violated on average. */
double rms_error = 0.0;
/* Sum of squared errors (in length units). */
double error_squared_sum = 0.0;
/* Largest squared error. */
double max_error_squared = 0.0;
} residual;
struct {
/* Time in seconds for the entire step. */
double step_total = 0.0;
/* Time in seconds to build the BVH tree of the surface.
* This is zero if the surface BVH was cached. */
double build_bvh = 0.0;
/* Time in seconds to find contact points, cumulative over substeps. */
double find_contacts = 0.0;
/* Time in seconds to solve constraints, cumulative over substeps. */
double solve_constraints = 0.0;
} timing;
};
private:
Params params_;
/** Length of each segment indexed by the index of the first point in the segment. */
Array<float> segment_lengths_cu_;
struct Contact {
float dist_sq_;
float3 normal_;
float3 point_;
};
Array<int> contacts_num_;
Array<Contact> contacts_;
/** Information about the most recent step solution. */
mutable Result result_;
public:
const Params &params() const;
Span<float> segment_lengths() const;
const Result &result() const;
void clear_result();
/* Initialize the solver for a given set of curves.
* The solver must be reinitialized if the curve set changes. */
void initialize(const Params &params,
const bke::CurvesGeometry &curves,
IndexMask curve_selection);
/* Solve constraints for an independent subset of curves. */
void step_curves(bke::CurvesGeometry &curves,
const Mesh *surface,
const bke::CurvesSurfaceTransforms &transforms,
Span<float3> start_positions,
IndexMask changed_curves,
bool update_error = false);
private:
void find_contact_points(const bke::CurvesGeometry &curves,
const Mesh *surface,
const bke::CurvesSurfaceTransforms &transforms,
float max_dist,
IndexMask changed_curves);
void apply_distance_constraint(float3 &point_a,
float3 &point_b,
float segment_length,
float weight_a,
float weight_b) const;
float get_distance_constraint_error(const float3 &point_a,
const float3 &point_b,
const float segment_length) const;
void apply_contact_constraint(float3 &point, float radius, const Contact &contact) const;
float get_contact_constraint_error(const float3 &point,
float radius,
const Contact &contact) const;
void solve_constraints(bke::CurvesGeometry &curves, IndexMask changed_curves) const;
void solve_curve_constraints(bke::CurvesGeometry &curves,
const VArray<float> &radius,
const IndexRange points) const;
void compute_error(const bke::CurvesGeometry &curves, IndexMask changed_curves) const;
void compute_curve_error(const bke::CurvesGeometry &curves,
const VArray<float> &radius,
const IndexRange points) const;
};
} // namespace blender::geometry

458
source/blender/geometry/intern/constraint_solver.cc
View File

@ -1,458 +0,0 @@
/* SPDX-License-Identifier: GPL-2.0-or-later */
/** \file
* \ingroup bke
*/
#include "GEO_constraint_solver.hh"
#include "BLI_index_mask_ops.hh"
#include "BKE_bvhutils.h"
#include "BKE_mesh.h"
#include "BKE_mesh_runtime.h"
#include "DNA_mesh_types.h"
#include "DNA_meshdata_types.h"
#include "DNA_object_types.h"
#include "PIL_time.h"
namespace blender::geometry {
using bke::CurvesGeometry;
using bke::CurvesSurfaceTransforms;
using threading::EnumerableThreadSpecific;
static const int curves_grain_size = 64;
const ConstraintSolver::Params &ConstraintSolver::params() const
{
return params_;
}
Span<float> ConstraintSolver::segment_lengths() const
{
return segment_lengths_cu_;
}
const ConstraintSolver::Result &ConstraintSolver::result() const
{
return result_;
}
void ConstraintSolver::clear_result()
{
result_ = Result{};
}
void ConstraintSolver::initialize(const Params &params,
const CurvesGeometry &curves,
IndexMask curve_selection)
{
params_ = params;
if (params_.use_length_constraints) {
segment_lengths_cu_.reinitialize(curves.points_num());
const Span<float3> positions_cu = curves.positions();
const OffsetIndices points_by_curve = curves.points_by_curve();
threading::parallel_for(curve_selection.index_range(), 256, [&](const IndexRange range) {
for (const int curve_i : curve_selection.slice(range)) {
const IndexRange points = points_by_curve[curve_i].drop_back(1);
float length_cu = 0.0f;
for (const int point_i : points) {
const float3 &p1_cu = positions_cu[point_i];
const float3 &p2_cu = positions_cu[point_i + 1];
length_cu = math::distance(p1_cu, p2_cu);
segment_lengths_cu_[point_i] = length_cu;
}
}
});
}
else {
segment_lengths_cu_.reinitialize(0);
}
if (params_.use_collision_constraints) {
contacts_num_.reinitialize(curves.points_num());
contacts_.reinitialize(params_.max_contacts_per_point * curves.points_num());
}
else {
contacts_num_.reinitialize(0);
contacts_.reinitialize(0);
}
}
void ConstraintSolver::step_curves(CurvesGeometry &curves,
const Mesh *surface,
const CurvesSurfaceTransforms &transforms,
const Span<float3> start_positions,
const IndexMask changed_curves,
bool update_error)
{
const double step_start = PIL_check_seconds_timer();
clear_result();
const OffsetIndices points_by_curve = curves.points_by_curve();
const MutableSpan<float3> positions = curves.positions_for_write();
const float max_substep_travel_distance = params_.max_travel_distance / params_.substep_count;
const float max_collision_distance = 2.0f * max_substep_travel_distance;
const float max_travel_distance_sq = params_.max_travel_distance * params_.max_travel_distance;
const bool clamp_travel = true;
/* Compute position delta per substep ("velocity") */
Array<float3> delta_substep(curves.points_num());
threading::parallel_for(
changed_curves.index_range(), curves_grain_size, [&](const IndexRange range) {
for (const int curve_i : changed_curves.slice(range)) {
const IndexRange points = points_by_curve[curve_i];
for (const int point_i : points) {
const float3 delta_step = positions[point_i] - start_positions[point_i];
positions[point_i] = start_positions[point_i];
const float travel_sq = len_squared_v3(delta_step);
if (travel_sq <= max_travel_distance_sq) {
delta_substep[point_i] = delta_step / params_.substep_count;
continue;
}
result_.max_travel_exceeded = true;
if (clamp_travel) {
delta_substep[point_i] = math::normalize(delta_step) * max_substep_travel_distance;
}
else {
delta_substep[point_i] = delta_step / params_.substep_count;
}
}
}
});
for ([[maybe_unused]] const int substep : IndexRange(params_.substep_count)) {
/* Set unconstrained position: x <- x + v*dt */
threading::parallel_for(
changed_curves.index_range(), curves_grain_size, [&](const IndexRange range) {
for (const int curve_i : changed_curves.slice(range)) {
const IndexRange points = points_by_curve[curve_i];
for (const int point_i : points) {
positions[point_i] += delta_substep[point_i];
}
}
});
if (params_.use_collision_constraints) {
find_contact_points(curves, surface, transforms, max_collision_distance, changed_curves);
}
solve_constraints(curves, changed_curves);
}
if (update_error) {
compute_error(curves, changed_curves);
}
result_.timing.step_total += PIL_check_seconds_timer() - step_start;
}
void ConstraintSolver::find_contact_points(const CurvesGeometry &curves,
const Mesh *surface,
const CurvesSurfaceTransforms &transforms,
const float max_dist,
const IndexMask changed_curves)
{
/* Should be set when initializing constraints */
BLI_assert(contacts_num_.size() == curves.points_num());
BLI_assert(contacts_.size() == curves.points_num() * params_.max_contacts_per_point);
contacts_num_.fill(0);
if (surface == nullptr) {
return;
}
const float curves_to_surface_scale = mat4_to_scale(transforms.curves_to_surface.ptr());
const float max_dist_su = curves_to_surface_scale * max_dist;
const float max_dist_sq_su = max_dist_su * max_dist_su;
const OffsetIndices points_by_curve = curves.points_by_curve();
const Span<MLoopTri> surface_looptris = surface->looptris();
const Span<float3> surface_positions = surface->vert_positions();
const Span<MLoop> surface_loops = surface->loops();
const double build_bvh_start = PIL_check_seconds_timer();
/** BVH tree of the surface mesh for finding collisions. */
BVHTreeFromMesh surface_bvh;
BKE_bvhtree_from_mesh_get(
&surface_bvh, surface, BVHTREE_FROM_LOOPTRI, params_.bvh_branching_factor);
BLI_SCOPED_DEFER([&]() { free_bvhtree_from_mesh(&surface_bvh); });
if (!surface_bvh.cached) {
result_.timing.build_bvh += PIL_check_seconds_timer() - build_bvh_start;
}
double find_contacts_start = PIL_check_seconds_timer();
threading::parallel_for(
changed_curves.index_range(), curves_grain_size, [&](const IndexRange range) {
for (const int curve_i : changed_curves.slice(range)) {
/* First point is anchored to the surface, ignore collisions. */
const IndexRange points = params_.use_root_constraints ?
points_by_curve[curve_i].drop_front(1) :
points_by_curve[curve_i];
for (const int point_i : points) {
const float3 &new_p = curves.positions()[point_i];
const MutableSpan<Contact> contacts = contacts_.as_mutable_span().slice(
params_.max_contacts_per_point * point_i, params_.max_contacts_per_point);
const float3 p_su = transforms.curves_to_surface * new_p;
BLI_bvhtree_range_query_cpp(
*surface_bvh.tree,
p_su,
max_dist_su,
[&](const int triangle_i, const float3 &co_su, float dist_sq_su) {
const MLoopTri &looptri = surface_looptris[triangle_i];
const float3 v0_su = surface_positions[surface_loops[looptri.tri[0]].v];
const float3 v1_su = surface_positions[surface_loops[looptri.tri[1]].v];
const float3 v2_su = surface_positions[surface_loops[looptri.tri[2]].v];
float3 closest_su;
closest_on_tri_to_point_v3(closest_su, co_su, v0_su, v1_su, v2_su);
dist_sq_su = len_squared_v3v3(co_su, closest_su);
if (dist_sq_su <= max_dist_sq_su) {
const int contacts_num = contacts_num_[point_i];
int insert_i;
if (contacts_num < params_.max_contacts_per_point) {
insert_i = contacts_num;
++contacts_num_[point_i];
}
else {
/* Replace the contact with the largest distance. */
insert_i = -1;
float max_dist_sq_su = dist_sq_su;
for (int contact_i : IndexRange(4)) {
if (contacts[contact_i].dist_sq_ > max_dist_sq_su) {
insert_i = contact_i;
max_dist_sq_su = contacts[contact_i].dist_sq_;
}
}
}
if (insert_i >= 0) {
float3 normal_su;
normal_tri_v3(normal_su, v0_su, v1_su, v2_su);
contacts[insert_i] = Contact{
max_dist_sq_su,
transforms.surface_to_curves_normal * float3{normal_su},
transforms.surface_to_curves * float3{closest_su}};
}
}
});
}
}
});
result_.timing.find_contacts += PIL_check_seconds_timer() - find_contacts_start;
}
void ConstraintSolver::apply_distance_constraint(float3 &point_a,
float3 &point_b,
const float segment_length,
const float weight_a,
const float weight_b) const
{
float length_cu;
const float3 direction = math::normalize_and_get_length(point_b - point_a, length_cu);
/* Constraint function */
const float C = length_cu - segment_length;
const float3 gradient = direction;
/* Lagrange multiplier */
const float lambda = -C / (weight_a + weight_b + params_.alpha);
point_a -= gradient * lambda * weight_a;
point_b += gradient * lambda * weight_b;
}
float ConstraintSolver::get_distance_constraint_error(const float3 &point_a,
const float3 &point_b,
const float segment_length) const
{
float length_cu = math::length(point_b - point_a);
/* Constraint function */
const float C = length_cu - segment_length;
return C;
}
void ConstraintSolver::apply_contact_constraint(float3 &point,
const float radius,
const ConstraintSolver::Contact &contact) const
{
/* Constraint function */
const float C = math::dot(point - contact.point_, contact.normal_) - radius;
const float3 gradient = contact.normal_;
/* Lagrange multiplier */
const float lambda = -C / (1.0f + params_.alpha);
if (lambda > 0.0f) {
point += gradient * lambda;
}
}
float ConstraintSolver::get_contact_constraint_error(
const float3 &point, const float radius, const ConstraintSolver::Contact &contact) const
{
/* Constraint function */
const float C = math::dot(point - contact.point_, contact.normal_) - radius;
return math::min(C, 0.0f);
}
void ConstraintSolver::solve_constraints(CurvesGeometry &curves,
const IndexMask changed_curves) const
{
const int solver_max_iterations = [&]() {
switch (params_.solver_type) {
case SolverType::Sequential:
return 1;
}
return 0;
}();
const double solve_constraints_start = PIL_check_seconds_timer();
const OffsetIndices points_by_curve = curves.points_by_curve();
VArray<float> radius = curves.attributes().lookup_or_default<float>(
"radius", ATTR_DOMAIN_POINT, 0.0f);
threading::parallel_for(
changed_curves.index_range(), curves_grain_size, [&](const IndexRange range) {
for (const int curve_i : changed_curves.slice(range)) {
const IndexRange points = points_by_curve[curve_i];
/* Solve constraints */
for ([[maybe_unused]] const int solver_i : IndexRange(solver_max_iterations)) {
solve_curve_constraints(curves, radius, points);
}
}
});
result_.timing.solve_constraints += PIL_check_seconds_timer() - solve_constraints_start;
}
void ConstraintSolver::solve_curve_constraints(CurvesGeometry &curves,
const VArray<float> &radius,
const IndexRange points) const
{
const int skipped_root_points = [&]() {
switch (params_.solver_type) {
/* The sequential solver only moves the 2nd point of each segment. */
case SolverType::Sequential:
return params_.use_root_constraints ? 1 : 0;
}
return 0;
}();
MutableSpan<float3> positions_cu = curves.positions_for_write();
result_.constraint_count = 0;
/* Distance constraints */
if (params_.use_length_constraints) {
result_.constraint_count += points.size() - 1;
for (const int segment_i : IndexRange(points.size() - 1)) {
const int point_a = points[segment_i];
const int point_b = points[segment_i + 1];
const float segment_length = segment_lengths_cu_[point_a];
float3 &pa = positions_cu[point_a];
float3 &pb = positions_cu[point_b];
const float w0 = (segment_i < skipped_root_points ? 0.0f : 1.0f);
const float w1 = 1.0f;
apply_distance_constraint(pa, pb, segment_length, w0, w1);
}
}
/* Contact constraints */
if (params_.use_collision_constraints) {
for (const int point_i : points.drop_front(skipped_root_points)) {
float3 &p = positions_cu[point_i];
const float radius_p = radius[point_i];
const int contacts_num = contacts_num_[point_i];
result_.constraint_count += contacts_num;
const Span<Contact> contacts = contacts_.as_span().slice(
params_.max_contacts_per_point * point_i, contacts_num);
for (const Contact &c : contacts) {
apply_contact_constraint(p, radius_p, c);
}
}
}
}
void ConstraintSolver::compute_error(const CurvesGeometry &curves,
const IndexMask changed_curves) const
{
const OffsetIndices points_by_curve = curves.points_by_curve();
VArray<float> radius = curves.attributes().lookup_or_default<float>(
"radius", ATTR_DOMAIN_POINT, 0.0f);
/* Accumulate error (no threading) */
for (const int curve_i : changed_curves) {
const IndexRange points = points_by_curve[curve_i];
compute_curve_error(curves, radius, points);
}
if (result_.constraint_count > 0) {
result_.residual.rms_error = sqrt(result_.residual.error_squared_sum /
result_.constraint_count);
}
else {
result_.residual.rms_error = 0.0;
}
if (result_.residual.max_error_squared > params_.error_threshold * params_.error_threshold) {
result_.status = Result::Status::ErrorNoConvergence;
}
}
void ConstraintSolver::compute_curve_error(const CurvesGeometry &curves,
const VArray<float> &radius,
const IndexRange points) const
{
Span<float3> positions_cu = curves.positions();
/* Distance constraints */
if (params_.use_length_constraints) {
for (const int segment_i : IndexRange(points.size() - 1)) {
const int point_a = points[segment_i];
const int point_b = points[segment_i + 1];
const float segment_length = segment_lengths_cu_[point_a];
const float3 &pa = positions_cu[point_a];
const float3 &pb = positions_cu[point_b];
const float error = get_distance_constraint_error(pa, pb, segment_length);
const double error_sq = error * error;
result_.residual.error_squared_sum += error_sq;
result_.residual.max_error_squared = std::max(result_.residual.max_error_squared, error_sq);
}
}
/* Contact constraints */
if (params_.use_collision_constraints) {
for (const int point_i : points) {
const float3 &p = positions_cu[point_i];
const float radius_p = radius[point_i];
const int contacts_num = contacts_num_[point_i];
const Span<Contact> contacts = contacts_.as_span().slice(
params_.max_contacts_per_point * point_i, contacts_num);
for (const Contact &c : contacts) {
const float error = get_contact_constraint_error(p, radius_p, c);
const double error_sq = error * error;
result_.residual.error_squared_sum += error_sq;
result_.residual.max_error_squared = std::max(result_.residual.max_error_squared,
error_sq);
}
}
}
}
} // namespace blender::geometry

86
source/blender/geometry/tests/GEO_constraint_solver_test.cc
View File

@ -1,86 +0,0 @@
// /* SPDX-License-Identifier: GPL-2.0-or-later */
// /** \file
// * \ingroup bke
// */
// #include "testing/testing.h"
// #include "BKE_curves.hh"
// #include "GEO_constraint_solver.hh"
// namespace blender::geometry::tests {
// using bke::CurvesGeometry;
// using bke::CurvesSurfaceTransforms;
// inline CurvesSurfaceTransforms create_curves_surface_transforms()
// {
// CurvesSurfaceTransforms transforms;
// transforms.curves_to_world = float4x4::identity();
// transforms.curves_to_surface = float4x4::identity();
// transforms.world_to_curves = float4x4::identity();
// transforms.world_to_surface = float4x4::identity();
// transforms.surface_to_world = float4x4::identity();
// transforms.surface_to_curves = float4x4::identity();
// transforms.surface_to_curves_normal = float4x4::identity();
// return transforms;
// }
// TEST(curves_constraints, LengthAndRootConstraint)
// {
// CurvesGeometry curves(9, 3);
// curves.fill_curve_types(CURVE_TYPE_POLY);
// const MutableSpan<int> point_offsets = curves.offsets_for_write();
// point_offsets[0] = 0;
// point_offsets[1] = 2;
// point_offsets[2] = 6;
// point_offsets[3] = 9;
// const MutableSpan<float3> positions = curves.positions_for_write();
// positions[0] = float3(.0f, 0, 0);
// positions[1] = float3(.1f, 0, 0);
// positions[2] = float3(.0f, 0, 1);
// positions[3] = float3(.1f, 0, 1);
// positions[4] = float3(.15f, 0, 1);
// positions[5] = float3(.0f, 0, 2);
// positions[6] = float3(.1f, 0, 2);
// positions[7] = float3(.15f, 0, 2);
// positions[8] = float3(.25f, 0, 2);
// const CurvesSurfaceTransforms transforms = create_curves_surface_transforms();
// ConstraintSolver solver;
// ConstraintSolver::Params params;
// params.use_collision_constraints = false;
// solver.initialize(params, curves, curves.curves_range());
// const Array<float3> orig_positions = curves.positions();
// positions[1].y += 0.1f;
// positions[3].y += 0.1f;
// positions[6].y += 0.1f;
// solver.step_curves(curves, nullptr, transforms, orig_positions, curves.curves_range());
// EXPECT_EQ(solver.result().status, ConstraintSolver::Result::Status::Ok);
// EXPECT_LE(solver.result().residual.max_error_squared, params.error_threshold *
// params.error_threshold); EXPECT_LE(solver.result().residual.rms_error,
// params.error_threshold);
// const float eps = 0.005f;
// EXPECT_V3_NEAR(positions[0], float3(0.0f, 0.0f, 0.0f), eps);
// EXPECT_V3_NEAR(positions[1], float3(0.065f, 0.075f, 0.0f), eps);
// EXPECT_V3_NEAR(positions[2], float3(0.0f, 0.0f, 1.0f), eps);
// EXPECT_V3_NEAR(positions[3], float3(0.0824425220, 0.057f, 1.0f), eps);
// EXPECT_V3_NEAR(positions[4], float3(0.118486673, 0.022f, 1.0f), eps);
// EXPECT_V3_NEAR(positions[5], float3(0.0f, 0.0f, 2.0f), eps);
// EXPECT_V3_NEAR(positions[6], float3(0.1f, 0.1f, 2.0f), eps);
// EXPECT_V3_NEAR(positions[7], float3(0.125f, 0.057f, 2.0f), eps);
// EXPECT_V3_NEAR(positions[8], float3(0.213f, 0.009f, 2.0f), eps);
// }
// } // namespace blender::geometry::tests

11
source/blender/geometry/tests/performance/CMakeLists.txt
View File

@ -1,11 +0,0 @@
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright 2014 Blender Foundation. All rights reserved.
set(INC
.
..
)
include_directories(${INC})
BLENDER_TEST_PERFORMANCE(GEO_constraint_solver_performance "bf_geometry")

364
source/blender/geometry/tests/performance/GEO_constraint_solver_performance_test.cc
View File

@ -1,364 +0,0 @@
/* SPDX-License-Identifier: GPL-2.0-or-later */
/** \file
* \ingroup bke
*/
// #include "testing/testing.h"
// #include "BLI_noise.hh"
// #include "BLI_rand.hh"
// #include "BKE_curves.hh"
// #include "BKE_idtype.h"
// #include "BKE_lib_id.h"
// #include "BKE_mesh.h"
// #include "DNA_mesh_types.h"
// #include "DNA_meshdata_types.h"
// #include "GEO_constraint_solver.hh"
// namespace blender::geometry::tests {
// using bke::CurvesGeometry;
// using bke::CurvesSurfaceTransforms;
// class CurveConstraintSolverPerfTestSuite : public testing::Test {
// public:
// /* Sets up Blender just enough to not crash on creating a mesh. */
// static void SetUpTestSuite()
// {
// testing::Test::SetUpTestSuite();
// //CLG_init();
// BLI_threadapi_init();
// //DNA_sdna_current_init();
// //BKE_blender_globals_init();
// BKE_idtype_init();
// }
// static void TearDownTestSuite()
// {
// BLI_threadapi_exit();
// //BKE_blender_atexit();
// //BKE_tempdir_session_purge();
// //BKE_appdir_exit();
// //CLG_exit();
// testing::Test::TearDownTestSuite();
// }
// };
// inline CurvesSurfaceTransforms create_curves_surface_transforms()
// {
// CurvesSurfaceTransforms transforms;
// transforms.curves_to_world = float4x4::identity();
// transforms.curves_to_surface = float4x4::identity();
// transforms.world_to_curves = float4x4::identity();
// transforms.world_to_surface = float4x4::identity();
// transforms.surface_to_world = float4x4::identity();
// transforms.surface_to_curves = float4x4::identity();
// transforms.surface_to_curves_normal = float4x4::identity();
// return transforms;
// }
// const uint32_t randomized_curves_seed = 12345;
// const uint32_t randomized_offset_seed = 23451;
// static float3 random_point_inside_unit_sphere(RandomNumberGenerator &rng)
// {
// const float theta = (float)(M_PI * 2.0) * rng.get_float();
// const float phi = acosf(2.0f * rng.get_float() - 1.0f);
// const float r = sqrt3f(rng.get_float());
// return r * float3{sinf(phi) * cosf(theta), sinf(phi) * sinf(theta), cosf(phi)};
// }
// static Mesh *create_noise_grid(const float noise = 0.1f, const int resolution = 100, const float
// size = 1.0f)
// {
// const int tot_verts = resolution * resolution;
// const int tot_polys = (resolution - 1) * (resolution - 1);
// const int tot_loops = tot_polys * 4;
// Mesh *mesh = BKE_mesh_new_nomain(tot_verts, 0, 0, tot_loops, tot_polys);
// const MutableSpan<MVert> mverts(mesh->mvert, mesh->totvert);
// const MutableSpan<MPoly> mpolys(mesh->mpoly, mesh->totpoly);
// const MutableSpan<MLoop> mloops(mesh->mloop, mesh->totloop);
// for (const int i : IndexRange(resolution)) {
// for (const int j : IndexRange(resolution)) {
// float3 co = float3((float)i - 0.5f, (float)j - 0.5f, 0) * size;
// const float3 offset = noise::perlin_float3_fractal_distorted(co, 2.0f, 0.5f, 0.0f);
// co += noise * (offset - float3(0.5f));
// const int vert_i = i * resolution + j;
// copy_v3_v3(mverts[vert_i].co, co);
// }
// }
// for (const int i : IndexRange(resolution - 1)) {
// for (const int j : IndexRange(resolution - 1)) {
// const int vert_i = i * resolution + j;
// const int poly_i = i * (resolution - 1) + j;
// const int loop_i = poly_i * 4;
// mpolys[poly_i].loopstart = loop_i;
// mpolys[poly_i].totloop = 4;
// mloops[loop_i + 0].v = vert_i;
// mloops[loop_i + 1].v = vert_i + 1;
// mloops[loop_i + 2].v = vert_i + resolution + 1;
// mloops[loop_i + 3].v = vert_i + resolution;
// }
// }
// BKE_mesh_calc_edges(mesh, false, false);
// return mesh;
// }
// static Mesh *create_torus(const int major_segments = 64,
// const int minor_segments = 16,
// const float major_radius = 1.0f,
// const float minor_radius = 0.5f)
// {
// const int tot_verts = major_segments * minor_segments;
// const int tot_polys = major_segments * minor_segments;
// const int tot_loops = tot_polys * 4;
// Mesh *mesh = BKE_mesh_new_nomain(tot_verts, 0, 0, tot_loops, tot_polys);
// const MutableSpan<MVert> mverts(mesh->mvert, mesh->totvert);
// const MutableSpan<MPoly> mpolys(mesh->mpoly, mesh->totpoly);
// const MutableSpan<MLoop> mloops(mesh->mloop, mesh->totloop);
// for (const int i : IndexRange(major_segments)) {
// for (const int j : IndexRange(minor_segments)) {
// const float alpha = 2.0 * M_PI * i / major_segments;
// const float beta = 2.0 * M_PI * j / minor_segments;
// const float3 tmp = float3(cosf(beta) + major_radius, 0.0f, sinf(beta)) * minor_radius;
// const float3 co = float3(tmp.x * cosf(alpha), tmp.x * sinf(alpha), tmp.z);
// const int vert_i = i * minor_segments + j;
// copy_v3_v3(mverts[vert_i].co, co);
// }
// }
// for (const int i : IndexRange(major_segments)) {
// for (const int j : IndexRange(minor_segments)) {
// const int vert_i = i * minor_segments + j;
// const int poly_i = i * minor_segments + j;
// const int loop_i = poly_i * 4;
// mpolys[poly_i].loopstart = loop_i;
// mpolys[poly_i].totloop = 4;
// const int di = i < major_segments - 1 ? minor_segments : minor_segments - tot_verts;
// const int dj = j < minor_segments - 1 ? 1 : 1 - minor_segments;
// mloops[loop_i + 0].v = vert_i;
// mloops[loop_i + 1].v = vert_i + dj;
// mloops[loop_i + 2].v = vert_i + di + dj;
// mloops[loop_i + 3].v = vert_i + di;
// }
// }
// BKE_mesh_calc_edges(mesh, false, false);
// return mesh;
// }
// static CurvesGeometry create_randomized_curves(const int curve_num,
// const int point_num_min,
// const int point_num_max,
// const float dist_min,
// const float dist_max)
// {
// RandomNumberGenerator rng(randomized_curves_seed);
// Array<int> point_offsets(curve_num + 1);
// const int point_num_range = std::max(point_num_max - point_num_min, 0);
// int point_num = 0;
// for (int curve_i : IndexRange(curve_num)) {
// point_offsets[curve_i] = point_num;
// point_num += point_num_min + (rng.get_int32() % point_num_range);
// }
// point_offsets.last() = point_num;
// CurvesGeometry curves(point_num, curve_num);
// curves.fill_curve_types(CURVE_TYPE_POLY);
// curves.offsets_for_write().copy_from(point_offsets);
// const MutableSpan<float3> positions = curves.positions_for_write();
// for (int curve_i : curves.curves_range()) {
// float3 pos = random_point_inside_unit_sphere(rng);
// for (int point_i : curves.points_for_curve(curve_i)) {
// positions[point_i] = pos;
// pos += rng.get_unit_float3() * interpf(dist_max, dist_min, rng.get_float());
// }
// }
// return curves;
// }
// static void randomized_point_offset(CurvesGeometry &curves,
// IndexMask changed_curves,
// float dist_min,
// float dist_max,
// float3 direction = float3(1, 0, 0),
// float cone_angle = M_PI)
// {
// RandomNumberGenerator rng(randomized_offset_seed);
// math::normalize(direction);
// float3 n1, n2;
// ortho_basis_v3v3_v3(n1, n2, direction);
// const MutableSpan<float3> positions = curves.positions_for_write();
// for (const int curve_i : changed_curves) {
// for (int point_i : curves.points_for_curve(curve_i)) {
// const float theta = rng.get_float() * cone_angle;
// const float phi = rng.get_float() * (float)(2.0 * M_PI);
// const float3 dir_p = direction * cosf(theta) + n1 * sinf(theta) * cosf(phi) +
// n2 * sinf(theta) * sinf(phi);
// const float dist_p = interpf(dist_max, dist_min, rng.get_float());
// positions[point_i] += dir_p * dist_p;
// }
// }
// }
// static void print_test_stats(const CurvesGeometry &curves, const ConstraintSolver::Result
// &result)
// {
// const testing::TestInfo *const test_info =
// testing::UnitTest::GetInstance()->current_test_info();
// std::cout << "Curves: " << curves.curves_num() << " Points: " << curves.points_num() <<
// std::endl; std::cout << "RMS=" << result.residual.rms_error
// << " max error=" << sqrt(result.residual.max_error_squared)
// << std::endl;
// std::cout << "total: " << result.timing.step_total * 1000.0f
// << " ms, build bvh: " << result.timing.build_bvh * 1000.0f
// << " ms, find contacts: " << result.timing.find_contacts * 1000.0f
// << " ms, solve: " << result.timing.solve_constraints * 1000.0f << " ms" <<
// std::endl;
// }
// template <typename ParamType>
// class SolverPerfTestSuite : public CurveConstraintSolverPerfTestSuite,
// public testing::WithParamInterface<ParamType> {
// public:
// struct TestResult {
// ParamType param;
// ConstraintSolver::Result solver_result;
// };
// static Vector<TestResult> &results()
// {
// static Vector<TestResult> results_;
// return results_;
// }
// static void SetUpTestSuite()
// {
// CurveConstraintSolverPerfTestSuite::SetUpTestSuite();
// results().clear();
// }
// static void TearDownTestSuite()
// {
// /* CSV printout for simple data import */
// std::cout << "Parameter,RMS Error,Max Error,Collision Detection (ms),Solve Time (ms)" <<
// std::endl; for (const TestResult &result : results()) {
// std::cout << result.param << "," << result.solver_result.residual.rms_error << ","
// << sqrtf(result.solver_result.residual.max_error_squared) << ","
// << (result.solver_result.timing.find_contacts * 1000.0) << ","
// << (result.solver_result.timing.solve_constraints * 1000.0) << std::endl;
// }
// CurveConstraintSolverPerfTestSuite::TearDownTestSuite();
// }
// void randomized_test(const ParamType &test_param,
// const ConstraintSolver::Params &solver_params,
// const int mesh_resolution = 100)
// {
// CurvesGeometry curves = create_randomized_curves(10000, 4, 50, 0.1f, 0.2f);
// const CurvesSurfaceTransforms transforms = create_curves_surface_transforms();
// Mesh *surface = create_torus(mesh_resolution * 2, mesh_resolution, 1.0f, 0.5f);
// ConstraintSolver solver;
// solver.initialize(solver_params, curves, curves.curves_range());
// const Array<float3> orig_positions = curves.positions();
// randomized_point_offset(curves, curves.curves_range(), 0.0f, 1.0f);
// solver.step_curves(curves, surface, transforms, orig_positions, curves.curves_range(),
// true);
// BKE_id_free(nullptr, surface);
// results().append(TestResult{test_param, solver.result()});
// }
// };
// using SolverIterationsTestSuite = SolverPerfTestSuite<int>;
// TEST_P(SolverIterationsTestSuite, RandomizedTest)
// {
// ConstraintSolver::Params params;
// params.solver_type = ConstraintSolver::SolverType::PositionBasedDynamics;
// params.max_solver_iterations = GetParam();
// randomized_test(GetParam(), params);
// }
// INSTANTIATE_TEST_SUITE_P(curves_constraints_performance,
// SolverIterationsTestSuite,
// testing::Values(1, 5, 10, 20, 50, 100),
// [](const testing::TestParamInfo<int> info) {
// std::stringstream ss;
// ss << "SolverIterationsTest_" << info.param;
// return ss.str();
// });
// using BVHBranchingFactorTestSuite = SolverPerfTestSuite<int>;
// TEST_P(BVHBranchingFactorTestSuite, RandomizedTest)
// {
// ConstraintSolver::Params params;
// params.bvh_branching_factor = GetParam();
// randomized_test(GetParam(), params);
// }
// INSTANTIATE_TEST_SUITE_P(curves_constraints_performance,
// BVHBranchingFactorTestSuite,
// testing::Values(2, 4, 6, 8, 16, 32),
// [](const testing::TestParamInfo<int> info) {
// std::stringstream ss;
// ss << "BVHBranchingFactorTest_" << info.param;
// return ss.str();
// });
// using MeshDensityTestSuite = SolverPerfTestSuite<int>;
// TEST_P(MeshDensityTestSuite, RandomizedTest)
// {
// ConstraintSolver::Params params;
// randomized_test(GetParam(), params, GetParam());
// }
// INSTANTIATE_TEST_SUITE_P(curves_constraints_performance,
// MeshDensityTestSuite,
// testing::Values(50, 100, 150, 200, 250),
// [](const testing::TestParamInfo<int> info) {
// std::stringstream ss;
// ss << "MeshDensityTest_" << info.param;
// return ss.str();
// });
// } // namespace blender::geometry::tests