Skip to content

Commit 4edb2f2

Browse files
HuangzizhouHuangzizhou
committed
Merge branch 'main' of https://github.com/polyfem/polyfem 
2 parents 9055589 + 0a76346 commit 4edb2f2

File tree

12 files changed

+236
-292
lines changed

12 files changed

+236
-292
lines changed

‎cmake/recipes/polyfem_data.cmake‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ else()
2828
PREFIX ${FETCHCONTENT_BASE_DIR}/polyfem-test-data
2929
SOURCE_DIR ${POLYFEM_DATA_DIR}
3030
GIT_REPOSITORY https://github.com/polyfem/polyfem-data
31-
GIT_TAG e90b99d13f030f3b4a59bd0997aa70fef66e8d73
31+
GIT_TAG 089ec98f9fdb9bf25fb413c13091102eb74685c5
3232
CONFIGURE_COMMAND ""
3333
BUILD_COMMAND ""
3434
INSTALL_COMMAND ""

‎cmake/recipes/polyfem_diff.cmake‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
# data (https://github.com/polyfem/polyfem-data)
2+
# License: MIT
3+
4+
if(TARGET polyfem::differentiable)
5+
return()
6+
endif()
7+
8+
include(ExternalProject)
9+
10+
set(POLYFEM_DIFF_DIR "${PROJECT_SOURCE_DIR}/diff-data/" CACHE PATH "Where should polyfem download and look for diff data?")
11+
12+
ExternalProject_Add(
13+
polyfem_diff_download
14+
PREFIX ${FETCHCONTENT_BASE_DIR}/polyfem-test-diff
15+
SOURCE_DIR ${POLYFEM_DIFF_DIR}
16+
GIT_REPOSITORY https://github.com/polyfem/differentiability-data
17+
GIT_TAG 7e49f248417987f0187a8bd0171954f486152a6d
18+
CONFIGURE_COMMAND ""
19+
BUILD_COMMAND ""
20+
INSTALL_COMMAND ""
21+
LOG_DOWNLOAD ON
22+
)
23+
24+
# Create a dummy target for convenience
25+
add_library(polyfem_diff INTERFACE)
26+
add_library(polyfem::diff ALIAS polyfem_diff)
27+
28+
add_dependencies(polyfem_diff polyfem_diff_download)
29+
30+
target_compile_definitions(polyfem_diff INTERFACE POLYFEM_DIFF_DIR=\"${POLYFEM_DIFF_DIR}\")

‎src/polyfem/assembler/FixedCorotational.cpp‎

Lines changed: 6 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -354,7 +354,6 @@ namespace polyfem::assembler
354354
}
355355

356356
const Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
357-
358357
const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
359358

360359
// Id + grad d
@@ -411,28 +410,23 @@ namespace polyfem::assembler
411410
}
412411

413412
Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
413+
const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
414414

415415
// Id + grad d
416-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
416+
def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
417417

418418
double lambda, mu;
419419
params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
420420

421421
Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp = compute_stiffness_from_def_grad(def_grad, lambda, mu);
422422

423-
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor(jac_it.size(), grad.size());
423+
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor = Eigen::Matrix<double, dim * dim, N>::Zero(jac_it.size(), grad.size());
424424

425425
for (size_t i = 0; i < local_disp.rows(); ++i)
426426
{
427-
for (size_t j = 0; j < local_disp.cols(); ++j)
428-
{
429-
Eigen::Matrix<double, dim, dim> temp(size(), size());
430-
temp.setZero();
431-
temp.row(j) = grad.row(i);
432-
temp = temp * jac_it;
433-
Eigen::Matrix<double, dim * dim, 1> temp_flattened(Eigen::Map<Eigen::Matrix<double, dim * dim, 1>>(temp.data(), temp.size()));
434-
delF_delU_tensor.col(i * size() + j) = temp_flattened;
435-
}
427+
for (size_t j = 0; j < dim; ++j)
428+
for (size_t k = 0; k < dim; ++k)
429+
delF_delU_tensor(dim * k + j, i * dim + j) = delF_delU(i, k);
436430
}
437431

438432
Eigen::Matrix<double, N, N> hessian = delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor;

‎src/polyfem/assembler/MooneyRivlin3ParamSymbolic.cpp‎

Lines changed: 6 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -437,9 +437,10 @@ namespace polyfem::assembler
437437
}
438438

439439
Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
440+
Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
440441

441442
// Id + grad d
442-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
443+
def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
443444
def_grad_T = def_grad.transpose();
444445

445446
const double t = 0;
@@ -451,7 +452,6 @@ namespace polyfem::assembler
451452
Eigen::Matrix<double, dim, dim> gradient_temp;
452453
autogen::generate_gradient_templated<dim>(c1, c2, c3, d1, def_grad_T, gradient_temp);
453454

454-
Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
455455
Eigen::Matrix<double, n_basis, dim> gradient = delF_delU * gradient_temp.transpose();
456456
G.noalias() += gradient * data.da(p);
457457
}
@@ -497,9 +497,10 @@ namespace polyfem::assembler
497497
}
498498

499499
Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
500+
const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
500501

501502
// Id + grad d
502-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
503+
def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
503504

504505
const double t = 0;
505506
const double c1 = c1_(data.vals.val.row(p), t, data.vals.element_id);
@@ -510,41 +511,15 @@ namespace polyfem::assembler
510511
Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp;
511512
autogen::generate_hessian_templated<dim>(c1, c2, c3, d1, def_grad, hessian_temp);
512513

513-
// Check by FD
514-
/*
515-
{
516-
double eps = 1e-7;
517-
Eigen::MatrixXd gradient_temp, gradient_temp_plus;
518-
autogen::generate_gradient(c1, c2, c3, d1, def_grad, gradient_temp);
519-
Eigen::MatrixXd x_ = def_grad;
520-
for (int j = 0; j < hessian_temp.cols(); ++j)
521-
{
522-
x_(j / dim, j % dim) += eps;
523-
autogen::generate_gradient(c1, c2, c3, d1, x_, gradient_temp_plus);
524-
Eigen::MatrixXd fd = (gradient_temp_plus - gradient_temp) / eps;
525-
526-
std::cout << "hess " << fd.transpose() << "\t" << hessian_temp.col(j).transpose() << std::endl;
527-
// for (int i = 0; i < hessian_temp.rows(); ++i)
528-
// {
529-
// double fd = (gradient_temp_plus(i) - gradient_temp(i)) / eps;
530-
// if (abs(fd - hessian_temp(i, j)) > 1e-7)
531-
// std::cout << "mismatch " << abs(fd - hessian_temp(i, j)) << std::endl;
532-
// }
533-
x_(j / dim, j % dim) -= eps;
534-
}
535-
}
536-
*/
537-
538-
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor(jac_it.size(), grad.size());
514+
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor = Eigen::Matrix<double, dim * dim, N>::Zero(jac_it.size(), grad.size());
539515

540516
for (size_t i = 0; i < local_disp.rows(); ++i)
541517
{
542518
for (size_t j = 0; j < local_disp.cols(); ++j)
543519
{
544520
Eigen::Matrix<double, dim, dim> temp(size(), size());
545521
temp.setZero();
546-
temp.row(j) = grad.row(i);
547-
temp = temp * jac_it;
522+
temp.row(j) = delF_delU.row(i);
548523
Eigen::Matrix<double, dim * dim, 1> temp_flattened(Eigen::Map<Eigen::Matrix<double, dim * dim, 1>>(temp.data(), temp.size()));
549524
delF_delU_tensor.col(i * size() + j) = temp_flattened;
550525
}

‎src/polyfem/assembler/NeoHookeanElasticity.cpp‎

Lines changed: 51 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -302,53 +302,85 @@ namespace polyfem::assembler
302302

303303
double NeoHookeanElasticity::compute_energy(const NonLinearAssemblerData &data) const
304304
{
305-
return compute_energy_aux<double>(data);
305+
if (size() == 2)
306+
{
307+
switch (data.vals.basis_values.size())
308+
{
309+
case 3:
310+
return compute_energy_aux<double, 3, 2>(data);
311+
case 6:
312+
return compute_energy_aux<double, 6, 2>(data);
313+
case 10:
314+
return compute_energy_aux<double, 10, 2>(data);
315+
default:
316+
return compute_energy_aux<double, Eigen::Dynamic, 2>(data);
317+
}
318+
}
319+
else // if (size() == 3)
320+
{
321+
assert(size() == 3);
322+
switch (data.vals.basis_values.size())
323+
{
324+
case 4:
325+
return compute_energy_aux<double, 4, 3>(data);
326+
case 10:
327+
return compute_energy_aux<double, 10, 3>(data);
328+
case 20:
329+
return compute_energy_aux<double, 20, 3>(data);
330+
default:
331+
return compute_energy_aux<double, Eigen::Dynamic, 3>(data);
332+
}
333+
}
306334
}
307335

308336
// Compute ∫ ½μ (tr(FᵀF) - 3 - 2ln(J)) + ½λ ln²(J) du
309-
template <typename T>
337+
template <typename T, int n_basis, int dim>
310338
T NeoHookeanElasticity::compute_energy_aux(const NonLinearAssemblerData &data) const
311339
{
312340
if constexpr (std::is_same_v<T, double>)
313341
{
314-
Eigen::MatrixXd local_disp(data.vals.basis_values.size(), size());
342+
Eigen::Matrix<T, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
315343
local_disp.setZero();
316344
for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
317345
{
318346
const auto &bs = data.vals.basis_values[i];
319347
for (size_t ii = 0; ii < bs.global.size(); ++ii)
320348
{
321-
for (int d = 0; d < size(); ++d)
349+
for (int d = 0; d < dim; ++d)
322350
{
323351
local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
324352
}
325353
}
326354
}
355+
327356
Eigen::VectorXd jacs;
328357
if (use_robust_jacobian)
329358
jacs = data.vals.eval_deformed_jacobian_determinant(data.x);
330-
Eigen::MatrixXd def_grad(size(), size());
359+
360+
Eigen::Matrix<T, dim, dim> def_grad(size(), size());
331361

332362
T energy = T(0.0);
333363

334364
const int n_pts = data.da.size();
335365
for (long p = 0; p < n_pts; ++p)
336366
{
337-
Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
367+
Eigen::Matrix<T, n_basis, dim> grad(data.vals.basis_values.size(), size());
338368

339369
for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
340370
{
341371
grad.row(i) = data.vals.basis_values[i].grad.row(p);
342372
}
343373

344-
const Eigen::MatrixXd jac_it = data.vals.jac_it[p];
374+
const Eigen::Matrix<T, dim, dim> jac_it = data.vals.jac_it[p];
345375
// Id + grad d
346-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::MatrixXd::Identity(size(), size());
376+
def_grad = (local_disp.transpose() * grad) * jac_it + Eigen::Matrix<T, dim, dim>::Identity(size(), size());
377+
347378
double lambda, mu;
348379
params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
380+
349381
const T J = use_robust_jacobian ? jacs(p) * jac_it.determinant() : def_grad.determinant();
350382
const T log_det_j = log(J);
351-
const T val = mu / 2 * ((def_grad.transpose() * def_grad).trace() - size() - 2 * log_det_j) +
383+
const T val = mu / 2 * (def_grad.squaredNorm() - size() - 2 * log_det_j) +
352384
lambda / 2 * log_det_j * log_det_j;
353385
energy += val * data.da(p);
354386
}
@@ -357,7 +389,7 @@ namespace polyfem::assembler
357389
else
358390
{
359391
typedef Eigen::Matrix<T, Eigen::Dynamic, 1> AutoDiffVect;
360-
typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
392+
typedef Eigen::Matrix<T, dim, dim> AutoDiffGradMat;
361393

362394
AutoDiffVect local_disp;
363395
get_local_disp(data, size(), local_disp);
@@ -372,7 +404,7 @@ namespace polyfem::assembler
372404
compute_disp_grad_at_quad(data, local_disp, p, size(), def_grad);
373405

374406
// Id + grad d
375-
for (int d = 0; d < size(); ++d)
407+
for (int d = 0; d < dim; ++d)
376408
def_grad(d, d) += T(1);
377409

378410
double lambda, mu;
@@ -432,7 +464,7 @@ namespace polyfem::assembler
432464
const auto &bs = data.vals.basis_values[i];
433465
for (size_t ii = 0; ii < bs.global.size(); ++ii)
434466
{
435-
for (int d = 0; d < size(); ++d)
467+
for (int d = 0; d < dim; ++d)
436468
{
437469
local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
438470
}
@@ -458,9 +490,10 @@ namespace polyfem::assembler
458490
}
459491

460492
Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
493+
Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
461494

462495
// Id + grad d
463-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
496+
def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
464497

465498
const double J = use_robust_jacobian ? jacs(p) * jac_it.determinant() : def_grad.determinant();
466499
const double log_det_j = log(J);
@@ -496,8 +529,6 @@ namespace polyfem::assembler
496529
double lambda, mu;
497530
params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
498531

499-
Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
500-
501532
Eigen::Matrix<double, dim, dim> gradient_temp = mu * def_grad - mu * (1 / J) * delJ_delF + lambda * log_det_j * (1 / J) * delJ_delF;
502533
Eigen::Matrix<double, n_basis, dim> gradient = delF_delU * gradient_temp.transpose();
503534

@@ -528,10 +559,7 @@ namespace polyfem::assembler
528559
const auto &bs = data.vals.basis_values[i];
529560
for (size_t ii = 0; ii < bs.global.size(); ++ii)
530561
{
531-
for (int d = 0; d < size(); ++d)
532-
{
533-
local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
534-
}
562+
local_disp.row(i) += bs.global[ii].val * data.x.block<dim, 1>(bs.global[ii].index * size(), 0).transpose();
535563
}
536564
}
537565

@@ -551,9 +579,10 @@ namespace polyfem::assembler
551579
}
552580

553581
Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
582+
const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
554583

555584
// Id + grad d
556-
def_grad = local_disp.transpose() * grad * jac_it + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
585+
def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
557586

558587
const double J = use_robust_jacobian ? jacs(p) * jac_it.determinant() : def_grad.determinant();
559588
double log_det_j = log(J);
@@ -606,16 +635,15 @@ namespace polyfem::assembler
606635

607636
Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp = (mu * id) + (((mu + lambda * (1 - log_det_j)) / (J * J)) * (g_j * g_j.transpose())) + (((lambda * log_det_j - mu) / (J)) * del2J_delF2);
608637

609-
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor(jac_it.size(), grad.size());
638+
Eigen::Matrix<double, dim * dim, N> delF_delU_tensor = Eigen::Matrix<double, dim * dim, N>::Zero(jac_it.size(), grad.size());
610639

611640
for (size_t i = 0; i < local_disp.rows(); ++i)
612641
{
613642
for (size_t j = 0; j < local_disp.cols(); ++j)
614643
{
615644
Eigen::Matrix<double, dim, dim> temp(size(), size());
616645
temp.setZero();
617-
temp.row(j) = grad.row(i);
618-
temp = temp * jac_it;
646+
temp.row(j) = delF_delU.row(i);
619647
Eigen::Matrix<double, dim * dim, 1> temp_flattened(Eigen::Map<Eigen::Matrix<double, dim * dim, 1>>(temp.data(), temp.size()));
620648
delF_delU_tensor.col(i * size() + j) = temp_flattened;
621649
}

‎src/polyfem/assembler/NeoHookeanElasticity.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@ namespace polyfem::assembler
7575
LameParameters params_;
7676

7777
// utility function that computes energy, the template is used for double, DScalar1, and DScalar2 in energy, gradient and hessian
78-
template <typename T>
78+
template <typename T, int n_basis, int dim>
7979
T compute_energy_aux(const NonLinearAssemblerData &data) const;
8080
template <int n_basis, int dim>
8181
void compute_energy_hessian_aux_fast(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const;

0 commit comments

Comments
 (0)