Integrate grad-grad operator.
780 {
782
784
788
791
792 auto &
locMat = AssemblyBoundaryEleOp::locMat;
793 locMat.resize(nb_rows, nb_cols,
false);
795
796 if (nb_cols && nb_rows) {
797
798 auto nb_gauss_pts = getGaussPts().size2();
799 auto t_w = getFTensor0IntegrationWeight();
800 auto t_coords = getFTensor1CoordsAtGaussPts();
801 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
802 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
803
804
805 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
808
810 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
811 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
812 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
813
814 auto next = [&]() {
815 ++t_w;
816 ++t_coords;
817 ++t_disp;
818 ++t_traction;
819 ++t_sdf_v;
820 ++t_grad_sdf_v;
821 ++t_hess_sdf_v;
822 ++t_normalized_normal;
823 };
824
825 auto face_data_vec_ptr =
827 auto face_gauss_pts_it = face_data_vec_ptr->begin();
828
830 auto nb_face_functions = row_data.
getN().size2() / 3;
832
833 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
834
835 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
836 face_data_vec_ptr);
837
838 auto check_face_contact = [&]() {
840 return false;
841
842 if (face_data_ptr) {
843 return true;
844 }
845 return false;
846 };
847
849
850#ifdef ENABLE_PYTHON_BINDING
852 if (ContactOps::sdfPythonWeakPtr.lock()) {
853 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
855 }
856#else
857 constexpr double c = 0;
858#endif
859
860 if (!
c && check_face_contact()) {
862 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
863 constexpr double eps = std::numeric_limits<float>::epsilon();
864 for (auto ii = 0; ii < 3; ++ii) {
866 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
867 t_spatial_coords_cx(ii) +=
eps * 1
i;
868 auto t_rhs_tmp =
871 for (int jj = 0; jj != 3; ++jj) {
872 auto v = t_rhs_tmp(jj).imag();
873 t_res_dU(jj, ii) =
v /
eps;
874 }
875 }
876
877 } else {
878
879#ifdef ENABLE_PYTHON_BINDING
880
881 if (ContactOps::sdfPythonWeakPtr.lock()) {
884
885 (-
c) * (t_hess_sdf_v(
i,
j) * t_grad_sdf_v(
k) * t_traction(
k) +
886 t_grad_sdf_v(
i) * t_hess_sdf_v(
k,
j) * t_traction(
k))
887
888 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
889
890 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
891 } else {
893 }
894#else
896#endif
897 }
898
899 auto alpha = t_w * getMeasure();
900
901 size_t rr = 0;
902 for (; rr != nb_rows / 3; ++rr) {
903
904 auto t_mat = getFTensor2FromArray<3, 3, 3>(
locMat, 3 * rr);
906
907 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
908 auto beta = alpha * t_row_base * t_col_base;
909 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
910 ++t_col_base;
911 ++t_mat;
912 }
913
914 ++t_row_base;
915 }
916 for (; rr < nb_face_functions; ++rr)
917 ++t_row_base;
918
919 next();
920 }
921 }
922
924}
#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)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
MatrixDouble locMat
local entity block matrix