727 {
729
731
735
736 auto nb_rows = row_data.getIndices().size();
737 auto nb_cols = col_data.getIndices().size();
738
739 auto &locMat = AssemblyBoundaryEleOp::locMat;
740 locMat.resize(nb_rows, nb_cols, false);
741 locMat.clear();
742
743 if (nb_cols && nb_rows) {
744
745 const size_t nb_gauss_pts = OP::getGaussPts().size2();
746
747 auto t_w = OP::getFTensor0IntegrationWeight();
748 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
749 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
750 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
751 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
752
753
754 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
757
758 auto t_sdf_v = getFTensor0FromVec(v_sdf);
759 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
760
761 auto next = [&]() {
762 ++t_w;
763 ++t_disp;
764 ++t_traction;
765 ++t_coords;
766 ++t_material_normal;
767 ++t_sdf_v;
768 ++t_grad_sdf_v;
769 };
770
771 auto face_data_vec_ptr =
773 auto face_gauss_pts_it = face_data_vec_ptr->begin();
774
775 auto t_row_base = row_data.getFTensor0N();
776 auto nb_face_functions = row_data.getN().size2();
777
779
780 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
781
782 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
783 face_data_vec_ptr);
784
785 auto check_face_contact = [&]() {
787 return false;
788
789 if (face_data_ptr) {
790 return true;
791 }
792 return false;
793 };
794
796
797#ifdef ENABLE_PYTHON_BINDING
799 if (ContactOps::sdfPythonWeakPtr.lock()) {
800 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
802 }
803#else
804 constexpr double c = 0;
805#endif
806
807 if (!
c && check_face_contact()) {
809 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
810 constexpr double eps = std::numeric_limits<float>::epsilon();
811 for (auto ii = 0; ii != 3; ++ii) {
813 t_traction(0), t_traction(1), t_traction(2)};
814 t_traction_cx(ii) +=
eps * 1
i;
815 auto t_rhs_tmp =
818 for (int jj = 0; jj != 3; ++jj) {
819 auto v = t_rhs_tmp(jj).imag();
820 t_res_dP(jj, ii) =
v /
eps;
821 }
822 }
823 } else {
824
825#ifdef ENABLE_PYTHON_BINDING
826 if (ContactOps::sdfPythonWeakPtr.lock()) {
828 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
830 t_res_dP(
i,
j) = t_cQ(
i,
j);
831 } else {
833 }
834#else
836#endif
837 }
838
839 const double alpha = t_w / 2.;
840 size_t rr = 0;
841 for (; rr != nb_rows / 3; ++rr) {
842
843 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
844 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
845
846 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
847 auto col_base = t_col_base(
i) * t_material_normal(
i);
848 const double beta = alpha * t_row_base * col_base;
849 t_mat(
i,
j) -= beta * t_res_dP(
i,
j);
850 ++t_col_base;
851 ++t_mat;
852 }
853
854 ++t_row_base;
855 }
856 for (; rr < nb_face_functions; ++rr)
857 ++t_row_base;
858
859 next();
860 }
861 }
862
864}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto checkSdf(EntityHandle fe_ent, std::map< int, Range > &sdf_map_range)
auto multiPointRhs(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 3 > &t_spatial_coords, FTensor::Tensor1< T3, 3 > &t_master_traction, MultiPointRhsType type, bool debug=false)
auto getSdf(OP_PTR op_ptr, MatrixDouble &contact_disp, int block_id, bool eval_hessian)
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
