521 {
523
528
529 const size_t nb_gauss_pts = getGaussPts().size2();
530
531#ifndef NDEBUG
534 "Wrong number of integration pts %ld != %ld",
536 }
537#endif
538
539 auto &nf = locF;
540 locF.clear();
541
542 auto t_w = getFTensor0IntegrationWeight();
543 auto t_coords = getFTensor1CoordsAtGaussPts();
544 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
545 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
546
547
548
549
550 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
553
555 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
556 auto t_normalize_normal = getFTensor1FromMat<3>(m_normals_at_pts);
557
558 auto next = [&]() {
559 ++t_w;
560 ++t_coords;
561 ++t_disp;
562 ++t_traction;
563 ++t_normalize_normal;
564 ++t_sdf_v;
565 ++t_grad_sdf_v;
566 };
567
568 auto face_data_vec_ptr =
570 auto face_gauss_pts_it = face_data_vec_ptr->begin();
571
572 auto nb_base_functions = data.getN().size2();
573 auto t_base = data.getFTensor0N();
574 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
575
577 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
578 face_data_vec_ptr);
579
580 auto check_face_contact = [&]() {
582 return false;
583
584 if (face_data_ptr) {
585 return true;
586 }
587 return false;
588 };
589
590#ifdef ENABLE_PYTHON_BINDING
592 if (ContactOps::sdfPythonWeakPtr.lock()) {
593 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
595 }
596#else
597 constexpr double c = 0;
598#endif
599
600 if (!
c && check_face_contact()) {
602 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
603 auto t_rhs_tmp =
multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
605 t_rhs(
i) = t_rhs_tmp(
i);
606
607 } else {
608
609#ifdef ENABLE_PYTHON_BINDING
611
612 if (ContactOps::sdfPythonWeakPtr.lock()) {
614 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
616 t_rhs(
i) = t_cQ(
i,
j) * t_traction(
j) +
617 (
c * inv_cn * t_sdf_v) * t_grad_sdf_v(
i);
618 } else {
619 t_rhs(
i) = t_traction(
i);
620 }
621#else
622 t_rhs(
i) = t_traction(
i);
623#endif
624 }
625
626 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
627 const double alpha = t_w * getMeasure();
628
629 size_t bb = 0;
630 for (; bb != nbRows / 3; ++bb) {
631 const double beta = alpha * t_base;
632 t_nf(
i) -= beta * t_rhs(
i);
633 ++t_nf;
634 ++t_base;
635 }
636 for (; bb < nb_base_functions; ++bb)
637 ++t_base;
638
639 next();
640 }
641
643}
#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.