521 auto nb_rows = row_data.getIndices().size();
522 auto nb_cols = col_data.getIndices().size();
525 locMat.resize(nb_rows, nb_cols,
false);
528 if (nb_cols && nb_rows) {
530 auto nb_gauss_pts = getGaussPts().size2();
531 auto t_w = getFTensor0IntegrationWeight();
532 auto t_coords = getFTensor1CoordsAtGaussPts();
533 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
534 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
537 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
542 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
543 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
544 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
554 ++t_normalized_normal;
557 auto face_data_vec_ptr =
559 auto face_gauss_pts_it = face_data_vec_ptr->begin();
561 auto t_row_base = row_data.getFTensor0N();
562 auto nb_face_functions = row_data.getN().size2() / 3;
565 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
567 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
570 auto check_face_contact = [&]() {
579 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
581 if (!
c && check_face_contact()) {
583 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
584 constexpr
double eps = std::numeric_limits<float>::epsilon();
585 for (
auto ii = 0; ii < 3; ++ii) {
587 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
588 t_spatial_coords_cx(ii) +=
eps * 1
i;
592 for (
int jj = 0; jj != 3; ++jj) {
593 auto v = t_rhs_tmp(jj).imag();
594 t_res_dU(jj, ii) =
v /
eps;
604 (-
c) * (t_hess_sdf_v(
i,
j) * t_grad_sdf_v(
k) * t_traction(
k) +
605 t_grad_sdf_v(
i) * t_hess_sdf_v(
k,
j) * t_traction(
k))
607 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
608 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
611 auto alpha = t_w * getMeasure();
614 for (; rr != nb_rows / 3; ++rr) {
616 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
617 auto t_col_base = col_data.getFTensor0N(gg, 0);
619 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
620 auto beta = alpha * t_row_base * t_col_base;
621 t_mat(
i,
j) += beta * t_res_dU(
i,
j);
628 for (; rr < nb_face_functions; ++rr)