M arap-mold.cpp => arap-mold.cpp +26 -41
@@ 23,8 23,7 @@ typedef Eigen::MatrixXi Derivedb;
typedef Eigen::MatrixXd Derivedbc;
typedef Eigen::MatrixXd DerivedU;
-Eigen::Quaternion<double> quat_angle_multiply(Eigen::Quaternion<double> qt,
- double f) {
+Eigen::Quaternion<double> quat_angle_multiply(Eigen::Quaternion<double> qt, double f) {
Eigen::AngleAxis<double> aa(qt);
aa.angle() *= f;
return Eigen::Quaternion<double>(aa);
@@ 32,10 31,8 @@ Eigen::Quaternion<double> quat_angle_multiply(Eigen::Quaternion<double> qt,
// Check whether `face` is moldable, given the plane
arap::Moldable::Enum arap::is_moldable(Eigen::Vector3d &remove_dir,
- std::vector<bool> &is_point_above,
- Eigen::MatrixXi &faces, int face,
- Eigen::Vector3d normal,
- double *out_dot) {
+ std::vector<bool> &is_point_above, Eigen::MatrixXi &faces,
+ int face, Eigen::Vector3d normal, double *out_dot) {
bool x_above = is_point_above[faces(face, 0)];
bool y_above = is_point_above[faces(face, 1)];
@@ 57,9 54,8 @@ arap::Moldable::Enum arap::is_moldable(Eigen::Vector3d &remove_dir,
return Moldable::Yes;
}
-bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
- const int dim, const Eigen::MatrixXi &b,
- arap::Data &data) {
+bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
+ const Eigen::MatrixXi &b, arap::Data &data) {
using namespace std;
using namespace Eigen;
typedef typename DerivedV::Scalar Scalar;
@@ 171,16 167,14 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
data.vel = MatrixXd::Zero(n, data.dim);
}
- return min_quad_with_fixed_precompute(Q, b, SparseMatrix<double>(), true,
- data.solver_data);
+ return min_quad_with_fixed_precompute(Q, b, SparseMatrix<double>(), true, data.solver_data);
}
// # Note: if ran with `use_dynamics` you must set `ARAPData::U0` to be your
// initial `U` before calling this. Otherwise it'll be 0, an you will have a bad
// time.
-void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
- arap::Data &data, int num_tets, Eigen::MatrixXi &TNeigh,
- Eigen::MatrixXd &U) {
+void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &data, int num_tets,
+ Eigen::MatrixXi &TNeigh, Eigen::MatrixXd &U) {
using namespace Eigen;
using namespace std;
assert(data.b.size() == bc.rows());
@@ 191,8 185,7 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
if (U.size() == 0) {
// terrible initial guess.. should at least copy input mesh
#ifndef NDEBUG
- cerr << "arap_solve: Using terrible initial guess for U. Try U = V."
- << endl;
+ cerr << "arap_solve: Using terrible initial guess for U. Try U = V." << endl;
#endif
U = MatrixXd::Zero(data.n, data.dim);
} else {
@@ 294,8 287,7 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
auto remove_dir = *data.mht_plane_n;
if (data.orientation->at(face_i) == Orientation::Below)
remove_dir = -remove_dir;
- Eigen::Vector3d rotation_axis =
- face_normal.cross(remove_dir).normalized();
+ Eigen::Vector3d rotation_axis = face_normal.cross(remove_dir).normalized();
auto x = remove_dir.dot(face_normal);
auto angle_between = atan2(1.0, x);
@@ 310,8 302,7 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
// multiple rotations.
// XXX(note): We're assuming this uses right-hand rotation, but the
// Eigen docs doens't say anything about it(!!).
- Eigen::Quaternion<double> qt(
- Eigen::AngleAxis<double>(angle_between, rotation_axis));
+ Eigen::Quaternion<double> qt(Eigen::AngleAxis<double>(angle_between, rotation_axis));
tet_rotations[tet_i].push_back(qt);
} break;
}
@@ 319,18 310,18 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
// Merge rotations for each tet
std::vector<Eigen::Quaternion<double>> merged_rotations(num_tets);
- for (int i = 0; i < num_tets; i++) {
+ for (int t = 0; t < num_tets; t++) {
Eigen::Quaternion<double> mold_quat;
- switch (tet_rotations[i].size()) {
+ switch (tet_rotations[t].size()) {
case 0: {
mold_quat.setIdentity();
} break;
case 1: {
- mold_quat = tet_rotations[i][0];
+ mold_quat = tet_rotations[t][0];
} break;
case 2: {
- Eigen::Quaternion<double> a(tet_rotations[i][0]);
- Eigen::Quaternion<double> b(tet_rotations[i][1]);
+ Eigen::Quaternion<double> a(tet_rotations[t][0]);
+ Eigen::Quaternion<double> b(tet_rotations[t][1]);
if (args.slerp_2_rots) {
mold_quat = a.slerp(0.5, b);
} else {
@@ 343,15 334,14 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
fprintf(stderr, "[WARN]: we're interpolating quats here (tet with >2 "
"boundary faces)\n");
Eigen::Quaternion<double> q(1, 0, 0, 0);
- double d = 1.0 / (double)tet_rotations[i].size();
- for (size_t ii = 0; ii < tet_rotations[i].size(); ii++) {
- tet_rotations[i][ii].w() *= d;
- q = tet_rotations[i][ii] * q;
+ double d = 1.0 / (double)tet_rotations[t].size();
+ for (size_t ii = 0; ii < tet_rotations[t].size(); ii++) {
+ q = quat_angle_multiply(tet_rotations[t][ii], d) * q;
}
mold_quat = q;
}
}
- merged_rotations[i] = mold_quat;
+ merged_rotations[t] = mold_quat;
}
// Blend rotations with neighboring tets
@@ 369,10 359,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
for (int n_i = 0; n_i < 4; n_i++) {
int neigh = TNeigh(tet_i, n_i);
if (neigh >= 0) {
- Eigen::Quaternion<double> neigh_qt =
- merged_rotations[TNeigh(tet_i, n_i)];
- neigh_qt = quat_angle_multiply(neigh_qt,
- args.blend / (double)num_neighbors);
+ Eigen::Quaternion<double> neigh_qt = merged_rotations[TNeigh(tet_i, n_i)];
+ neigh_qt = quat_angle_multiply(neigh_qt, args.blend / (double)num_neighbors);
qt = neigh_qt * qt;
}
}
@@ 399,16 387,14 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
} else {
eff_R.resize(Rdim, num_rots * Rdim);
for (int r = 0; r < num_rots; r++) {
- eff_R.block(0, Rdim * r, Rdim, Rdim) =
- Rotations.block(0, Rdim * data.G(r), Rdim, Rdim);
+ eff_R.block(0, Rdim * r, Rdim, Rdim) = Rotations.block(0, Rdim * data.G(r), Rdim, Rdim);
}
}
MatrixXd Dl;
if (data.with_dynamics) {
- assert(
- data.M.rows() == n &&
- "No mass matrix. Call arap_precomputation if changing with_dynamics");
+ assert(data.M.rows() == n &&
+ "No mass matrix. Call arap_precomputation if changing with_dynamics");
const double h = data.h;
assert(h != 0);
// Dl = 1./(h*h*h)*M*(-2.*V0 + Vm1) - fext;
@@ 435,8 421,7 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
if (bc.size() > 0) {
bcc = bc.col(c);
}
- bool ok =
- igl::min_quad_with_fixed_solve(data.solver_data, Bc, bcc, Beq, Uc);
+ bool ok = igl::min_quad_with_fixed_solve(data.solver_data, Bc, bcc, Beq, Uc);
if (!ok) {
fprintf(stderr, "`igl::min_quad_with_fixed_solve` returned err!!\n");
}
M arap-mold.h => arap-mold.h +10 -13
@@ 88,8 88,8 @@ struct Data {
Eigen::MatrixXd U0;
Data()
- : n(0), G(), energy(igl::ARAP_ENERGY_TYPE_DEFAULT), with_dynamics(false),
- f_ext(), h(1), ym(1), max_iter(10), K(), CSM(), solver_data(), b(),
+ : n(0), G(), energy(igl::ARAP_ENERGY_TYPE_DEFAULT), with_dynamics(false), f_ext(), h(1),
+ ym(1), max_iter(10), K(), CSM(), solver_data(), b(),
dim(-1) // force this to be set by _precomputation
{};
};
@@ 100,9 100,7 @@ struct Args {
bool slerp_2_rots; /* SLERP rotations to join 2. Otherwise multiply by 0.5 */
bool rotate_before; /* Rotate normal before checking moldability */
- Args()
- : should_blend(false), blend(0.0), slerp_2_rots(false),
- rotate_before(false){};
+ Args() : should_blend(false), blend(0.0), slerp_2_rots(false), rotate_before(false){};
};
// Compute necessary information to start using an ARAP deformation
@@ 115,19 113,18 @@ struct Args {
// b #b list of "boundary" fixed vertex indices into V
// Outputs:
// data struct containing necessary precomputation
-bool precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
- const int dim, const Eigen::MatrixXi &b, Data &data);
+bool precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
+ const Eigen::MatrixXi &b, Data &data);
// Inputs:
// bc #b by dim list of boundary conditions
// data struct containing necessary precomputation and parameters
// U #V by dim initial guess
-void iter(const Args &, const Eigen::MatrixXd &bc, Data &data, int num_tets,
- Eigen::MatrixXi &, Eigen::MatrixXd &U);
+void iter(const Args &, const Eigen::MatrixXd &bc, Data &data, int num_tets, Eigen::MatrixXi &,
+ Eigen::MatrixXd &U);
-Moldable::Enum is_moldable(Eigen::Vector3d &remove_dir,
- std::vector<bool> &is_point_above,
- Eigen::MatrixXi &faces, int face,
- Eigen::Vector3d normal, double *out_dot);
+Moldable::Enum is_moldable(Eigen::Vector3d &remove_dir, std::vector<bool> &is_point_above,
+ Eigen::MatrixXi &faces, int face, Eigen::Vector3d normal,
+ double *out_dot);
} // namespace arap
M libigl => libigl +1 -1
@@ 1,1 1,1 @@
-Subproject commit aed32a7322b9047576149c483b058bfde36c1eb1
+Subproject commit f6b406427400ed7ddb56cfc2577b6af571827c8c
M main.cpp => main.cpp +64 -192
@@ 31,7 31,6 @@ void ensure_moldability();
void update_rest_state();
void remesh();
int ray_intersect(Vector3d, Vector3d, MatrixXd &, MatrixXi &);
-void compute_forces(MatrixXd &, MatrixXi &);
void compute_things(MatrixXd &, MatrixXi &, bool = false);
int intersect3D_RayTriangle(Vector3d, Vector3d, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
@@ 49,6 48,7 @@ MatrixXi TF2T; // Face to tets mapping
MatrixXi TT2F; // Tet to faces mapping
MatrixXi TNeigh; // Neighbor tets for each tet
+MatrixXd RestState;
// Indices of the fixed vertices in TV.
VectorXi boundary(BOUNDARY_SIZE);
// ARAP'd vertices
@@ 69,13 69,14 @@ std::vector<double> dots;
arap::Data arap_data;
arap::Args arap_args;
-
// Whether to stitch together unmoldable geometry
bool stitching = false;
// Toggle running of ARAP
bool run_arap = false;
// Show rotations for all *faces*
bool show_list = false;
+// Visualize the forces on each vertex when stitching
+bool show_forces = false;
// Cut plane and point with `float`, for ImGui.
Eigen::Vector3f fplane_p;
@@ 98,11 99,13 @@ void update_vectors() {
}
void update_rest_state() {
- arap::precomputation(U, TT, TV.cols(), boundary, arap_data);
+ RestState = U;
+ arap::precomputation(RestState, TT, RestState.cols(), boundary, arap_data);
}
// Draw the custom menu
void draw_menu() {
+ bool recompute = false;
bool px = ImGui::SliderFloat("plane_p.x", &fplane_p[0], mins.x(), maxs.x(), "%.4f");
bool py = ImGui::SliderFloat("plane_p.y", &fplane_p[1], mins.y(), maxs.y(), "%.4f");
bool pz = ImGui::SliderFloat("plane_p.z", &fplane_p[2], mins.z(), maxs.z(), "%.4f");
@@ 115,8 118,7 @@ void draw_menu() {
plane_p = fplane_p.cast<double>();
plane_n = fplane_n.cast<double>().normalized();
update_vectors();
- draw_plane(viewer);
- compute_things(U, TF, true);
+ recompute |= true;
}
if (!run_arap) {
@@ 128,29 130,32 @@ void draw_menu() {
}
if (ImGui::Button("Reset")) {
- U = TV;
- compute_things(U, TF);
+ U = RestState;
+ recompute |= true;
}
if (ImGui::Button("Update Rest State")) {
update_rest_state();
+ recompute |= true;
}
if (ImGui::Button("Step ARAP")) {
step_arap();
}
- if (ImGui::Button("Compute Forces")) {
- compute_forces(U, TF);
- }
+ ImGui::Checkbox("Show Forces", &show_forces);
ImGui::Checkbox("Blend rotations", &arap_args.should_blend);
if (arap_args.should_blend) {
ImGui::SliderFloat("blend", &arap_args.blend, 0.0, 1.0, "%.2f");
}
+ if (ImGui::Checkbox("Stiching", &stitching)) {
+ if (!stitching)
+ Forces = MatrixXd::Zero(U.rows(), 3);
+ recompute |= true;
+ }
- ImGui::Checkbox("Stiching", &stitching);
ImGui::Checkbox("SLERP 2 rotations", &arap_args.slerp_2_rots);
ImGui::Checkbox("Rotate before mold check", &arap_args.rotate_before);
@@ 183,6 188,9 @@ void draw_menu() {
}
ImGui::EndChild();
}
+ if (recompute) {
+ compute_things(U, TF);
+ }
}
// Handle key down events.
@@ 192,7 200,7 @@ bool key_down(Viewer &viewer, unsigned char key, int mod) {
} else if (key == '.') {
step_arap();
} else if (key == 'R') {
- U = TV;
+ RestState = U = TV;
compute_things(U, TF);
}
return false;
@@ 202,15 210,10 @@ bool key_down(Viewer &viewer, unsigned char key, int mod) {
void step_arap() {
Eigen::MatrixXd bc(BOUNDARY_SIZE, 3);
for (int i = 0; i < BOUNDARY_SIZE; i++)
- bc.row(i) = TV.row(boundary(i));
+ bc.row(i) = RestState.row(boundary(i));
arap::iter(arap_args, bc, arap_data, TT.rows(), TNeigh, U);
compute_things(U, TF);
- if (stitching) {
- compute_forces(U, TF);
- } else {
- Forces = MatrixXd::Zero(U.rows(), 3);
- }
}
// Compute all the things that we need each time the vertex set changes. Here's
@@ 275,19 278,51 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
}
viewer.data().set_colors(colors);
}
-}
-// Here we look at the faces that are not moldable, and try to match them up with nearby faces; by
-// joining up faces like this we hope to fill in cavities, such that the final object is moldalbe.
-void merge_step() {}
+ if (stitching) {
+ Forces = MatrixXd::Zero(vertices.rows(), 3);
+ std::vector<bool> seen(vertices.rows());
+ for (int f = 0; f < faces.rows(); f++) {
+ auto m = moldable[f];
+ // XXX(todo): what do we do about `Crossing`? For now, skip as well.
+ if (m != arap::Moldable::No)
+ continue;
+ Eigen::Vector3d ray_dir =
+ (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
+
+ for (int i = 0; i < 3; i++) {
+ int v = faces(f, i);
+ // Only do a vertex one time; this is okay since the direction is +- plane_n, and not the
+ // triangle normal.
+ if (seen[v])
+ continue;
+ seen[v] = true;
+
+ Vector3d vertex = vertices.row(v);
+
+ int other_i = ray_intersect(ray_dir, vertex, vertices, faces);
+ if (other_i == -1)
+ continue;
+
+ Vector3d other_vertex = vertices.row(other_i);
+ auto force_direction = other_vertex - vertex;
+ // XXX(todo): what should the magnitude of the force be? Should probably take a closer look
+ // at the dynamics implementation in `arap-mold.cpp`.
+ (void)force_direction; /* unused warning */
+ Forces.row(v) = force_direction;
+ }
+ }
+ arap_data.f_ext = Forces;
+ }
+}
// This is called before each draw call. Compute normals and set right colors.
// XXX(note): much of this should probably just be done after `step_arap`, since
// the rest is pretty constant.
bool pre_draw(Viewer &viewer) {
+ draw_plane(viewer);
if (run_arap) {
step_arap();
- compute_things(U, TF);
}
return false;
}
@@ 296,8 331,7 @@ bool pre_draw(Viewer &viewer) {
void draw_plane(Viewer &viewer) {
// Must clear the previous plane from the screen. Maybe there's a better way
// of doing this than to clear the entire thing.
- viewer.data().clear();
- viewer.data().set_mesh(TV, TF);
+ viewer.data().lines.resize(0, 0);
Eigen::Vector3d dir1(plane_n.cross(Eigen::Vector3d(0, 0, 1)));
if (dir1.norm() < 0.5) {
@@ 348,173 382,12 @@ void draw_plane(Viewer &viewer) {
Eigen::RowVector3d c(1.0, 1.0, 1.0);
viewer.data().add_edges(planeEA, planeEB, c);
-}
-// Deform the object so that moldability is ensured.
-void ensure_moldability() {
- // We'll do this by finding all faces that are not moldable, and change their
- // positions in `U` such that their faces would be moldable. This is weird,
- // because altering the position of a vertex which is adjacent to a
- // non-moldable face may make another face adjacent to that vertex, which was
- // previsously moldable, non-moldable. It'll be something like this:
- //
- // q = make_queue()>
- // for face in faces {
- // if non-moldable(face) {
- // q.insert(face)
- // }
- // }
- // while q is not empty {
- // f <- pop(q)
- // v <- offending_vertex(f)
- // Rt = mold_rotation(f, removal_dir)
- // normal = Rt * normal(f);
- // p = argmin(p)(v.dot(normal))
- // U[v] += d * normal; /* This is not right, but something like this */
- // }
- //
- Eigen::MatrixXd normals;
- // XXX(todo): move this to after all places in which U is modified. We pretty
- // much always need to have this up to date.
- igl::per_face_normals(U, TF, normals);
-
- std::vector<bool> is_point_above(U.rows());
- auto above_plane = [&](Eigen::Vector3d v) { return (v - plane_p).dot(plane_n) >= 0.0; };
- for (int i = 0; i < U.rows(); i++) {
- is_point_above[i] = above_plane(U.row(i));
- }
-
- auto furthest_vertex = [](Eigen::Vector3d dir, int f) {
- double x = dir.dot(U.row(TF(f, 0)));
- double y = dir.dot(U.row(TF(f, 1)));
- double z = dir.dot(U.row(TF(f, 2)));
- if (y <= x && z <= x)
- return 0;
- if (x <= y && z <= y)
- return 1;
- if (x <= z && y <= z)
- return 2;
-
- { /* XXX(debug) */
- std::cerr << dir << std::endl << std::endl;
- std::cerr << U.row(TF(f, 0)) << std::endl << std::endl;
- std::cerr << U.row(TF(f, 1)) << std::endl << std::endl;
- std::cerr << U.row(TF(f, 2)) << std::endl << std::endl;
- fprintf(stderr, "%lf %lf %lf\n", x, y, z);
- }
- assert(false && "This is impossible");
- return -1;
- };
-
- auto closest_vertex = [](Eigen::Vector3d dir, int f) {
- double x = dir.dot(U.row(TF(f, 0)));
- double y = dir.dot(U.row(TF(f, 1)));
- double z = dir.dot(U.row(TF(f, 2)));
- if (x <= y && x <= z)
- return 0;
- if (y <= x && y <= z)
- return 1;
- if (z <= x && z <= y)
- return 2;
- { /* XXX(debug) */
- std::cerr << dir << std::endl << std::endl;
- std::cerr << U.row(TF(f, 0)) << std::endl << std::endl;
- std::cerr << U.row(TF(f, 1)) << std::endl << std::endl;
- std::cerr << U.row(TF(f, 2)) << std::endl << std::endl;
- fprintf(stderr, "%lf %lf %lf\n", x, y, z);
- }
- assert(false && "This is impossible (closest_vertex)");
- return -1;
- };
-
- // auto highest_vertex = [](int f) {
- // double x = plane_n.dot(U.row(TF(f, 0)));
- // double y = plane_n.dot(U.row(TF(f, 1)));
- // double z = plane_n.dot(U.row(TF(f, 2)));
- // if (x >= y && x >= z)
- // return 0;
- // if (y >= x && y >= z)
- // return 1;
- // if (z >= x && z >= y)
- // return 2;
- // };
-
- std::vector<int> indices(TF.rows());
- for (int i = 0; i < TF.rows(); i++) {
- auto normal = normals.row(i);
- Eigen::Vector3d removal_dir(plane_n);
- auto moldable = arap::is_moldable(removal_dir, is_point_above, TF, i, normal, NULL);
- if (moldable == arap::Moldable::No) {
- indices.push_back(i);
- }
- }
-
- while (indices.size() > 0) {
- // We have to drag all points towards the point that's highest along
- // `plane_n` in order to make it moldable.
- int face = indices.back();
- indices.pop_back();
- Eigen::Vector3d normal = normals.row(face);
-
- Eigen::Vector3d removal_dir(plane_n);
- auto moldable = arap::is_moldable(removal_dir, is_point_above, TF, face, normal, NULL);
- // XXX(todo): how to handle crossing faces?
- if (moldable != arap::Moldable::No)
- continue;
- auto drag_dir = (normal - plane_n.dot(normal) * plane_n).normalized();
- { /* XXX(debug): NaN check */
- if (drag_dir[0] != drag_dir[0]) {
- std::cerr << drag_dir << std::endl;
- std::cerr << normal << std::endl;
- std::cerr << plane_n << std::endl;
- std::cerr << plane_n.dot(normal) << std::endl;
- std::cerr << plane_n.dot(normal) * plane_n << std::endl;
- }
- }
-
- auto closest = closest_vertex(removal_dir, face);
- U.row(TF(face, closest)) += drag_dir * 0.01;
-
- // auto top = highest_vertex(face);
- auto far = furthest_vertex(drag_dir, face);
- (void)far;
-
- // XXX(todo): check neighboring faces of these vertices and add to `indices`
- // as needed.
+ if (show_forces) {
+ viewer.data().add_edges(U, U + Forces, Eigen::RowVector3d(1.0, 0.0, 1.0));
}
}
-void compute_forces(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces) {
- Forces.resize(vertices.rows(), 3);
- for (int f = 0; f < TF.rows(); f++) {
- auto m = moldable[f];
- // XXX(todo): what do we do about `Crossing`? For now, skip as well.
- if (m != arap::Moldable::No)
- continue;
- Eigen::Vector3d ray_dir = (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
-
- for (int i = 0; i < 3; i++) {
- int v_i = faces(f, i);
- // XXX(todo): this will check each vertex once per non-moldable adjacent face, when we only
- // need to check once.
- Vector3d vertex = vertices.row(v_i);
-
- int other_i = ray_intersect(ray_dir, vertex, vertices, faces);
- if (other_i == -1)
- continue;
-
- Vector3d other_vertex = vertices.row(other_i);
- auto force_direction = other_vertex - vertex;
- // XXX(todo): what should the magnitude of the force be? Should probably
- // take a closer look at the dynamics implementation in `arap-mold.cpp`.
- (void)force_direction; /* unused warning */
- // XXX(todo): these forces can be visualized with libigl (?)
- Forces.row(v_i) = force_direction;
- }
- }
- arap_data.f_ext = Forces;
-}
-
// Shoot a ray from `p` in direction `v`, and see if it intersects with any of
// the `faces`. If it does, return the index of the closest vertex to the
// intersection point. If not, return `-1`.
@@ 640,7 513,8 @@ int main(int argc, char *argv[]) {
TF2T = tetmesh_out.face2tets;
TT2F = tetmesh_out.tet2faces;
TNeigh = tetmesh_out.neighbors;
- U = TV;
+ RestState = U = TV;
+ Forces = MatrixXd::Zero(TV.rows(), 3);
boundary << 0, 1, 2, 3; // 19, 22, 121, 122, 221, 224, 321, 326;
@@ 669,7 543,6 @@ int main(int argc, char *argv[]) {
arap_data.h = 0.01;
arap::precomputation(TV, TT, TV.cols(), boundary, arap_data);
-
maxs = V.colwise().maxCoeff().cast<float>();
mins = V.colwise().minCoeff().cast<float>();
scale = (maxs - mins).minCoeff();
@@ 689,9 562,8 @@ int main(int argc, char *argv[]) {
viewer.data().clear();
viewer.data().set_mesh(TV, TF);
viewer.data().set_face_based(true);
- viewer.core.is_animating = true;
+ viewer.core().is_animating = true;
- draw_plane(viewer);
compute_things(TV, TF);
viewer.launch();
M mesh.cpp => mesh.cpp +2 -2
@@ 57,8 57,7 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
}
::tetrahedralize(&tb, &tetgen_in, &tetgen_out);
} catch (int e) {
- cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__
- << ": " << e << endl;
+ cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__ << ": " << e << endl;
return 1;
}
@@ 66,6 65,7 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
if (tetgen_out.numberoftetrahedra == 0) {
cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__
<< ": failed to create any tets!" << endl;
+ return 1;
}
// Vertices