662 auto nb_rows = row_data.getIndices().size();
663 auto nb_cols = col_data.getIndices().size();
666 locMat.resize(nb_rows, nb_cols,
false);
669 if (nb_cols && nb_rows) {
671 const size_t nb_gauss_pts = getGaussPts().size2();
673 auto t_w = getFTensor0IntegrationWeight();
674 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
675 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
676 auto t_coords = getFTensor1CoordsAtGaussPts();
677 auto t_material_normal = getFTensor1NormalsAtGaussPts();
680 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
685 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
697 auto face_data_vec_ptr =
699 auto face_gauss_pts_it = face_data_vec_ptr->begin();
701 auto t_row_base = row_data.getFTensor0N();
702 auto nb_face_functions = row_data.getN().size2();
706 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
708 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
711 auto check_face_contact = [&]() {
720 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
722 if (!
c && check_face_contact()) {
724 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
725 constexpr
double eps = std::numeric_limits<float>::epsilon();
726 for (
auto ii = 0; ii != 3; ++ii) {
728 t_traction(0), t_traction(1), t_traction(2)};
729 t_traction_cx(ii) +=
eps * 1
i;
733 for (
int jj = 0; jj != 3; ++jj) {
734 auto v = t_rhs_tmp(jj).imag();
735 t_res_dP(jj, ii) =
v /
eps;
741 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
743 t_res_dP(
i,
j) = t_cQ(
i,
j);
746 const double alpha = t_w / 2.;
748 for (; rr != nb_rows / 3; ++rr) {
750 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
751 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
753 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
754 auto col_base = t_col_base(
i) * t_material_normal(
i);
755 const double beta = alpha * t_row_base * col_base;
756 t_mat(
i,
j) += beta * t_res_dP(
i,
j);
763 for (; rr < nb_face_functions; ++rr)