Class dedicated to integrate operator.
528 {
530
535
536 const size_t nb_gauss_pts = getGaussPts().size2();
537
538#ifndef NDEBUG
541 "Wrong number of integration pts %ld != %ld",
543 }
544#endif
545
548
549 auto t_w = getFTensor0IntegrationWeight();
550 auto t_coords = getFTensor1CoordsAtGaussPts();
551 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
552 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
553
554
555
556
557 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
560
562 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
563 auto t_normalize_normal = getFTensor1FromMat<3>(m_normals_at_pts);
564
565 auto next = [&]() {
566 ++t_w;
567 ++t_coords;
568 ++t_disp;
569 ++t_traction;
570 ++t_normalize_normal;
571 ++t_sdf_v;
572 ++t_grad_sdf_v;
573 };
574
575 auto face_data_vec_ptr =
577 auto face_gauss_pts_it = face_data_vec_ptr->begin();
578
579 auto nb_base_functions = data.
getN().size2();
581 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
582
584 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
585 face_data_vec_ptr);
586
587 auto check_face_contact = [&]() {
589 return false;
590
591 if (face_data_ptr) {
592 return true;
593 }
594 return false;
595 };
596
597#ifdef ENABLE_PYTHON_BINDING
599 if (ContactOps::sdfPythonWeakPtr.lock()) {
600 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
602 }
603#else
604 constexpr double c = 0;
605#endif
606
607 if (!
c && check_face_contact()) {
609 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
610 auto t_rhs_tmp =
multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
612 t_rhs(
i) = t_rhs_tmp(
i);
613
614 } else {
615
616#ifdef ENABLE_PYTHON_BINDING
618
619 if (ContactOps::sdfPythonWeakPtr.lock()) {
621 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
623 t_rhs(
i) = t_cQ(
i,
j) * t_traction(
j) +
624 (
c * inv_cn * t_sdf_v) * t_grad_sdf_v(
i);
625 } else {
626 t_rhs(
i) = t_traction(
i);
627 }
628#else
629 t_rhs(
i) = t_traction(
i);
630#endif
631 }
632
633 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
634 const double alpha = t_w * getMeasure();
635
636 size_t bb = 0;
637 for (; bb !=
nbRows / 3; ++bb) {
638 const double beta = alpha * t_base;
639 t_nf(
i) -= beta * t_rhs(
i);
640 ++t_nf;
641 ++t_base;
642 }
643 for (; bb < nb_base_functions; ++bb)
644 ++t_base;
645
646 next();
647 }
648
650}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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)
FTensor::Index< 'l', 3 > l
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.
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....
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows