712 auto nb_rows = row_data.getIndices().size();
713 auto nb_cols = col_data.getIndices().size();
716 locMat.resize(nb_rows, nb_cols,
false);
719 if (nb_cols && nb_rows) {
721 const size_t nb_gauss_pts = OP::getGaussPts().size2();
723 auto t_w = OP::getFTensor0IntegrationWeight();
724 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
725 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
726 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
727 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
730 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
735 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
747 auto face_data_vec_ptr =
749 auto face_gauss_pts_it = face_data_vec_ptr->begin();
751 auto t_row_base = row_data.getFTensor0N();
752 auto nb_face_functions = row_data.getN().size2();
756 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
758 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
761 auto check_face_contact = [&]() {
771 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
774 constexpr
double c = 0;
777 if (!
c && check_face_contact()) {
779 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
780 constexpr
double eps = std::numeric_limits<float>::epsilon();
781 for (
auto ii = 0; ii != 3; ++ii) {
783 t_traction(0), t_traction(1), t_traction(2)};
784 t_traction_cx(ii) +=
eps * 1
i;
788 for (
int jj = 0; jj != 3; ++jj) {
789 auto v = t_rhs_tmp(jj).imag();
790 t_res_dP(jj, ii) =
v /
eps;
797 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
799 t_res_dP(
i,
j) = t_cQ(
i,
j);
805 const double alpha = t_w / 2.;
807 for (; rr != nb_rows / 3; ++rr) {
809 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
810 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
812 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
813 auto col_base = t_col_base(
i) * t_material_normal(
i);
814 const double beta = alpha * t_row_base * col_base;
815 t_mat(
i,
j) += beta * t_res_dP(
i,
j);
822 for (; rr < nb_face_functions; ++rr)