555 auto nb_rows = row_data.getIndices().size();
556 auto nb_cols = col_data.getIndices().size();
558 auto &locMat = AssemblyBoundaryEleOp::locMat;
559 locMat.resize(nb_rows, nb_cols,
false);
562 if (nb_cols && nb_rows) {
564 auto nb_gauss_pts = getGaussPts().size2();
565 auto t_w = getFTensor0IntegrationWeight();
566 auto t_coords = getFTensor1CoordsAtGaussPts();
567 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
568 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
571 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
576 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
577 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
578 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
588 ++t_normalized_normal;
591 auto face_data_vec_ptr =
593 auto face_gauss_pts_it = face_data_vec_ptr->begin();
595 auto t_row_base = row_data.getFTensor0N();
596 auto nb_face_functions = row_data.getN().size2() / 3;
599 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
601 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
604 auto check_face_contact = [&]() {
614 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
617 constexpr
double c = 0;
620 if (!
c && check_face_contact()) {
622 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
623 constexpr
double eps = std::numeric_limits<float>::epsilon();
624 for (
auto ii = 0; ii < 3; ++ii) {
626 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
627 t_spatial_coords_cx(ii) +=
eps * 1
i;
631 for (
int jj = 0; jj != 3; ++jj) {
632 auto v = t_rhs_tmp(jj).imag();
633 t_res_dU(jj, ii) =
v /
eps;
645 (-
c) * (t_hess_sdf_v(
i,
j) * t_grad_sdf_v(
k) * t_traction(
k) +
646 t_grad_sdf_v(
i) * t_hess_sdf_v(
k,
j) * t_traction(
k))
648 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
650 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
657 auto alpha = t_w * getMeasure();
660 for (; rr != nb_rows / 3; ++rr) {
662 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
663 auto t_col_base = col_data.getFTensor0N(gg, 0);
665 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
666 auto beta = alpha * t_row_base * t_col_base;
667 t_mat(
i,
j) += beta * t_res_dU(
i,
j);
674 for (; rr < nb_face_functions; ++rr)