Skip to content

prototype of a ESP improvement/fix #281

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
fsichetti wants to merge 1 commit into from
Closed

prototype of a ESP improvement/fix #281

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions python/src/barrier/barrier.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,13 @@ void define_barrier(py::module_& m)
)ipc_Qu8mg5v7")
.def(py::init());

// PROTOTYPE: needed to test the per-barrier mollifier power from Python.
py::class_<
InversePowerBarrier<>, Barrier, std::shared_ptr<InversePowerBarrier<>>>(
m, "InversePowerBarrier", "Inverse power barrier h(d) / d^p.")
.def(py::init<double>(), py::arg("power"))
.def_property_readonly("power", &InversePowerBarrier<>::power);

py::class_<CubicBarrier<>, Barrier, std::shared_ptr<CubicBarrier<>>>(
m, "CubicBarrier", R"ipc_Qu8mg5v7(
Cubic barrier function from [Ando 2024].
Expand Down
1 change: 1 addition & 0 deletions python/src/potentials/barrier_potential.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,6 +277,7 @@ void define_esp_potential(py::module& m)
.def_readonly("quad_order", &ESPParameters::quad_order)
.def_readonly("area_weights", &ESPParameters::area_weights)
.def_readwrite("local_quadrature_transfer", &ESPParameters::local_quadrature_transfer)
.def_readwrite("barrier", &ESPParameters::barrier) // PROTOTYPE
.def_readonly("integration_type", &ESPParameters::integration_type);

py::class_<ESPPotential>(m, "ESPPotential")
Expand Down
27 changes: 22 additions & 5 deletions src/ipc/esp/esp_local_transfer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,15 @@ namespace {
struct Element {
std::vector<index_t> vertices, edges;
double weight;
double ee_scale = 1.0; // PROTOTYPE: area-consistent EE weight factor
};
// PROTOTYPE: reference length for area-consistent EE weights (0 = off).
double area_ref(const ESPParameters& p)
{
const auto& cfg = ee_mollifier_config();
return cfg.area_ref_dbar ? p.dbar : cfg.area_ref;
}

Element element(const CollisionMesh& mesh, const ESPParameters& p, size_t k)
{
Element out;
Expand All @@ -88,11 +96,15 @@ namespace {
out.edges.push_back(mesh.faces_to_edges()(k, i));
}
out.weight = p.area_weights ? mesh.face_areas()(k) / 9.0 : 1.0;
out.ee_scale = ee_area_scale_factor(
mesh.face_areas()(k) / 3.0, area_ref(p), /*face=*/true);
} else if ((k -= mesh.num_faces()) < mesh.num_codim_edges()) {
const auto e = mesh.codim_edges()[k];
out.vertices = { mesh.edges()(e, 0), mesh.edges()(e, 1) };
out.edges = { e };
out.weight = p.area_weights ? mesh.edge_area(e) / 4.0 : 1.0;
out.ee_scale = ee_area_scale_factor(
mesh.edge_area(e) / 2.0, area_ref(p), /*face=*/false);
} else {
const auto v = mesh.codim_vertices()[k - mesh.num_codim_edges()];
out.vertices = { v };
Expand Down Expand Up @@ -231,12 +243,16 @@ namespace {
if (const auto* b = dynamic_cast<const InversePowerBarrier<>*>(
params.barrier.get()))
power = std::max(1, int(std::lround(b->power())) + 1);
if (ee_mollifier_config().order_override > 0) // PROTOTYPE
power = ee_mollifier_config().order_override;

for (const auto& pair : pairs) {
const std::array<index_t, 4> prim = { mesh.edges()(pair.ea, 0),
mesh.edges()(pair.ea, 1),
mesh.edges()(pair.eb, 0),
mesh.edges()(pair.eb, 1) };
const std::array<index_t, 4> prim = { {
mesh.edges()(pair.ea, 0),
mesh.edges()(pair.ea, 1),
mesh.edges()(pair.eb, 0),
mesh.edges()(pair.eb, 1),
} };
Eigen::Matrix<T, 4, 3> positions;
if constexpr (Order == 0) {
for (int i = 0; i < 4; ++i)
Expand All @@ -257,14 +273,15 @@ namespace {
X.row(prim[0]), X.row(prim[1]), X.row(prim[2]), X.row(prim[3]),
value(d2));
T m =
Math<T>::cubic_spline(sqrt(d2) / params.ee_support()) * 1.5
Math<T>::cubic_spline(sqrt(d2) / params.ee_support()) * ee_mollifier_config().ee_weight /* PROTOTYPE */
* edge_edge_mollifier<T>(
positions.row(0).transpose(), positions.row(1).transpose(),
positions.row(2).transpose(), positions.row(3).transpose(),
types, d2);
const T initial = m;
for (int k = 1; k < power; ++k)
m = m * initial;
m = m * elem.ee_scale; // PROTOTYPE
const Eigen::Vector3<T> q =
(positions.row(0) + u * (positions.row(1) - positions.row(0)))
.transpose();
Expand Down
4 changes: 3 additions & 1 deletion src/ipc/esp/esp_parameters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#pragma once

#include <ipc/barrier/barrier.hpp>
#include <ipc/gcp/distance/mollifier_config.hpp>
#include <ipc/utils/logger.hpp>

#include <algorithm>
Expand Down Expand Up @@ -95,7 +96,8 @@ struct ESPParameters {
/// Requires use_near_far=true and the default vertex quadrature rule.
/// Currently supports frictionless contact only.
/// Has no effect in 2D or when dbar_factor=0. Disabled by default.
bool local_quadrature_transfer = false;
bool local_quadrature_transfer =
ee_mollifier_config().local_transfer_default; // PROTOTYPE env default
const IntegrationType integration_type;

/// @brief Near/far split parameter dbar / dhat, in [0, 2].
Expand Down
35 changes: 32 additions & 3 deletions src/ipc/esp/esp_potential.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ namespace {
// - InversePowerBarrier(p): order = round(p) + 1, so p=1 -> 2, p=2 -> 3.
int mollifier_order_for_barrier(const std::shared_ptr<Barrier>& barrier)
{
if (ee_mollifier_config().order_override > 0) { // PROTOTYPE
return ee_mollifier_config().order_override;
}
if (const auto* ip =
dynamic_cast<const InversePowerBarrier<>*>(barrier.get())) {
const int p = static_cast<int>(std::lround(ip->power()));
Expand Down Expand Up @@ -125,6 +128,26 @@ namespace {
elem.weight = params.area_weights ? mesh.vertex_area(v) : 1.;
return elem;
}

// PROTOTYPE: make an EE sample stand for about one r x r patch (length r
// on rods) instead of a fixed fraction of its element. Rest geometry only.
double ee_area_scale(
const CollisionMesh& mesh,
const IntegrationElement& elem,
const ESPParameters& params)
{
const auto& cfg = ee_mollifier_config();
const double r = cfg.area_ref_dbar ? params.dbar : cfg.area_ref;
if (elem.face >= 0) {
return ee_area_scale_factor(
mesh.face_areas()(elem.face) / 3.0, r, /*face=*/true);
}
if (elem.num_edges == 1) {
return ee_area_scale_factor(
mesh.edge_area(elem.edges[0]) / 2.0, r, /*face=*/false);
}
return 1.0;
}
} // namespace

double ESPPotential::operator()(
Expand Down Expand Up @@ -275,7 +298,7 @@ double ESPPotential::operator()(
double mollifier =
Math<double>::cubic_spline(
dist / m_params.ee_support())
* 1.5;
* ee_mollifier_config().ee_weight /* PROTOTYPE */;
mollifier *= edge_edge_mollifier<double>(
X.row(ea).transpose(),
X.row(eb).transpose(),
Expand All @@ -285,6 +308,8 @@ double ESPPotential::operator()(
mollifier,
mollifier_order_for_barrier(
m_params.barrier));
mollifier = mollifier
* ee_area_scale(mesh, elem, m_params); // PROTOTYPE

if (use_nf) {
const double p_near = PointPotentialHelper::
Expand Down Expand Up @@ -593,7 +618,7 @@ Eigen::VectorXd ESPPotential::gradient(

T mollifier = Math<T>::cubic_spline(
dist / m_params.ee_support())
* 1.5;
* ee_mollifier_config().ee_weight /* PROTOTYPE */;
mollifier *= edge_edge_mollifier<T>(
positionsT.row(0).transpose(),
positionsT.row(1).transpose(),
Expand All @@ -604,6 +629,8 @@ Eigen::VectorXd ESPPotential::gradient(
mollifier,
mollifier_order_for_barrier(
m_params.barrier));
mollifier = mollifier
* ee_area_scale(mesh, elem, m_params); // PROTOTYPE

const ESPCollisionDict<PointType::EDGE>& dict =
*(iter->second);
Expand Down Expand Up @@ -1022,7 +1049,7 @@ Eigen::SparseMatrix<double> ESPPotential::hessian(

T mollifier = Math<T>::cubic_spline(
dist / m_params.ee_support())
* 1.5;
* ee_mollifier_config().ee_weight /* PROTOTYPE */;
mollifier *= edge_edge_mollifier<T>(
positionsT.row(0).transpose(),
positionsT.row(1).transpose(),
Expand All @@ -1033,6 +1060,8 @@ Eigen::SparseMatrix<double> ESPPotential::hessian(
mollifier,
mollifier_order_for_barrier(
m_params.barrier));
mollifier = mollifier
* ee_area_scale(mesh, elem, m_params); // PROTOTYPE

const ESPCollisionDict<PointType::EDGE>& dict =
*(iter->second);
Expand Down
4 changes: 2 additions & 2 deletions src/ipc/esp/esp_tangential_collisions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,7 @@ void TangentialCollisions::build(
vertices.row(e11).transpose(), dist_sqr);
double mol =
Math<double>::cubic_spline(dist / params.ee_support())
* 1.5;
* ee_mollifier_config().ee_weight /* PROTOTYPE */;
mol *= edge_edge_mollifier<double>(
vertices.row(e00).transpose(),
vertices.row(e01).transpose(),
Expand Down Expand Up @@ -461,7 +461,7 @@ void TangentialCollisions::build(
vertices.row(e10).transpose(), vertices.row(e11).transpose(),
dist_sqr_ee);
double mollifier =
Math<double>::cubic_spline(dist_ee / params.ee_support()) * 1.5;
Math<double>::cubic_spline(dist_ee / params.ee_support()) * ee_mollifier_config().ee_weight /* PROTOTYPE */;
mollifier *= edge_edge_mollifier<double>(
vertices.row(e00).transpose(), vertices.row(e01).transpose(),
vertices.row(e10).transpose(), vertices.row(e11).transpose(),
Expand Down
32 changes: 32 additions & 0 deletions src/ipc/gcp/distance/mollifier.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,6 +296,38 @@ std::array<HeavisideType, 4> edge_edge_mollifier_type(
const double dist_sqr)
{
std::array<HeavisideType, 4> mtypes {};
if (ee_mollifier_config().smoothstep) {
const double ell2 = ee_mollifier_config().rel_width
? dist_sqr * MOLLIFIER_THRESHOLD_EPS
: ee_mollifier_config().sqrt_width
? std::sqrt(dist_sqr) * MOLLIFIER_THRESHOLD_EPS
* ee_mollifier_config().L
: sqr(ee_mollifier_config().ell);
const auto t = [&](const double pe) {
return pe - dist_sqr >= ell2 ? HeavisideType::ONE
: HeavisideType::VARIANT;
};
mtypes[0] = t(point_edge_distance(ea0, eb0, eb1));
mtypes[1] = t(point_edge_distance(ea1, eb0, eb1));
mtypes[2] = t(point_edge_distance(eb0, ea0, ea1));
mtypes[3] = t(point_edge_distance(eb1, ea0, ea1));
return mtypes;
}
if (ee_mollifier_config().fixed || ee_mollifier_config().sqrt_width) {
const double denom = ee_mollifier_config().fixed
? sqr(ee_mollifier_config().ell)
: std::sqrt(dist_sqr) * MOLLIFIER_THRESHOLD_EPS
* ee_mollifier_config().L;
const auto t = [&](const double pe) {
return pe - dist_sqr >= denom ? HeavisideType::ONE
: HeavisideType::VARIANT;
};
mtypes[0] = t(point_edge_distance(ea0, eb0, eb1));
mtypes[1] = t(point_edge_distance(ea1, eb0, eb1));
mtypes[2] = t(point_edge_distance(eb0, ea0, ea1));
mtypes[3] = t(point_edge_distance(eb1, ea0, ea1));
return mtypes;
}
mtypes[0] = (point_edge_distance(ea0, eb0, eb1) - dist_sqr)
>= dist_sqr * MOLLIFIER_THRESHOLD_EPS
? HeavisideType::ONE
Expand Down
45 changes: 43 additions & 2 deletions src/ipc/gcp/distance/mollifier.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "mollifier.hpp"

#include <ipc/distance/point_line.hpp>
#include <ipc/gcp/distance/mollifier_config.hpp>

namespace ipc {
template <typename scalar, int dim>
Expand All @@ -28,8 +29,48 @@ scalar edge_edge_mollifier(
const std::array<HeavisideType, 4>& mtypes,
const scalar& dist_sqr)
{
const scalar da = dist_sqr * MOLLIFIER_THRESHOLD_EPS;
const scalar db = dist_sqr * MOLLIFIER_THRESHOLD_EPS;
if (ee_mollifier_config().smoothstep) {
// PROTOTYPE: H(c / l) with H(t) = 10t^3 - 15t^4 + 6t^5 (C2), where
// c = sqrt(d_PE^2 - d^2) is the endpoint clearance and l is fixed.
using std::sqrt;
const scalar ell2 = ee_mollifier_config().rel_width
? scalar(dist_sqr * MOLLIFIER_THRESHOLD_EPS)
: ee_mollifier_config().sqrt_width
? scalar(sqrt(dist_sqr)
* (MOLLIFIER_THRESHOLD_EPS * ee_mollifier_config().L))
: scalar(ee_mollifier_config().ell * ee_mollifier_config().ell);
const auto H = [&](const scalar& c2) -> scalar {
if (c2 <= 0) {
return scalar(0.);
}
if (c2 >= ell2) {
return scalar(1.);
}
const scalar t = sqrt(c2 / ell2);
return t * t * t * (10. + t * (-15. + 6. * t));
};
using PE = PointEdgeDistance<scalar, 3>;
scalar w(1.);
if (mtypes[0] == HeavisideType::VARIANT)
w = w * H(PE::point_edge_sqr_distance(ea0, eb0, eb1) - dist_sqr);
if (mtypes[1] == HeavisideType::VARIANT)
w = w * H(PE::point_edge_sqr_distance(ea1, eb0, eb1) - dist_sqr);
if (mtypes[2] == HeavisideType::VARIANT)
w = w * H(PE::point_edge_sqr_distance(eb0, ea0, ea1) - dist_sqr);
if (mtypes[3] == HeavisideType::VARIANT)
w = w * H(PE::point_edge_sqr_distance(eb1, ea0, ea1) - dist_sqr);
return w;
}
scalar denom = dist_sqr * MOLLIFIER_THRESHOLD_EPS;
if (ee_mollifier_config().fixed) { // PROTOTYPE
denom = scalar(ee_mollifier_config().ell * ee_mollifier_config().ell);
} else if (ee_mollifier_config().sqrt_width) { // PROTOTYPE
using std::sqrt;
denom = sqrt(dist_sqr)
* (MOLLIFIER_THRESHOLD_EPS * ee_mollifier_config().L);
}
const scalar da = denom;
const scalar db = denom;
scalar a = (mtypes[0] == HeavisideType::VARIANT)
? Math<scalar>::mollifier(
(PointEdgeDistance<scalar, 3>::point_edge_sqr_distance(
Expand Down
Loading
Loading