8#ifndef __CONTACTOPS_HPP__
9#define __CONTACTOPS_HPP__
14struct CommonData :
public boost::enable_shared_from_this<CommonData> {
19 static SmartPetscObj<Vec>
23 constexpr int ghosts[] = {0, 1, 2};
25 createGhostVector(m_field.
get_comm(),
48 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
53 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
64 SDFPython() =
default;
65 virtual ~SDFPython() =
default;
67 MoFEMErrorCode sdfInit(
const std::string py_file) {
72 auto main_module = bp::import(
"__main__");
73 mainNamespace = main_module.attr(
"__dict__");
74 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
76 sdfFun = mainNamespace[
"sdf"];
77 sdfGradFun = mainNamespace[
"grad_sdf"];
78 sdfHessFun = mainNamespace[
"hess_sdf"];
80 }
catch (bp::error_already_set
const &) {
90 py_list_to_std_vector(
const boost::python::object &iterable) {
91 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
92 boost::python::stl_input_iterator<T>());
97 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
105 sdf = bp::extract<double>(sdfFun(
t, x, y, z, tx, ty, tz));
107 }
catch (bp::error_already_set
const &) {
117 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
118 std::vector<double> &grad_sdf
126 py_list_to_std_vector<double>(sdfGradFun(
t, x, y, z, tx, ty, tz));
128 }
catch (bp::error_already_set
const &) {
143 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
144 std::vector<double> &hess_sdf
152 py_list_to_std_vector<double>(sdfHessFun(
t, x, y, z, tx, ty, tz));
154 }
catch (bp::error_already_set
const &) {
168 bp::object mainNamespace;
170 bp::object sdfGradFun;
171 bp::object sdfHessFun;
174static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
180 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz)>;
183 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz)>;
186 boost::function<FTensor::Tensor2_symmetric<double, 3>(
187 double t,
double x,
double y,
double z,
double tx,
double ty,
191 double tx,
double ty,
double tz) {
194 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
197 "Failed python call");
206 double tx,
double ty,
double tz) {
208 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
209 std::vector<double> grad_sdf;
211 "Failed python call");
220 double tx,
double ty,
double tz) {
222 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
223 std::vector<double> hess_sdf;
225 "Failed python call");
227 hess_sdf[2], hess_sdf[3],
228 hess_sdf[4], hess_sdf[5]};
234template <
int DIM, IntegrationType I,
typename BoundaryEleOp>
237template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
240template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
243template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
246template <
int DIM,
typename BoundaryEleOp>
250 boost::shared_ptr<CommonData> common_data_ptr);
257template <
int DIM,
typename AssemblyBoundaryEleOp>
261 boost::shared_ptr<CommonData> common_data_ptr);
262 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
272template <
int DIM,
typename AssemblyBoundaryEleOp>
276 const std::string col_field_name,
277 boost::shared_ptr<CommonData> common_data_ptr);
278 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
279 EntitiesFieldData::EntData &col_data);
290template <
int DIM,
typename AssemblyBoundaryEleOp>
294 const std::string row_field_name,
const std::string col_field_name,
295 boost::shared_ptr<CommonData> common_data_ptr);
296 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
297 EntitiesFieldData::EntData &col_data);
309 template <
int DIM, IntegrationType I>
318 template <
int DIM, IntegrationType I>
322 template <
int DIM, IntegrationType I>
326 template <
int DIM, IntegrationType I>
341inline double w(
const double sdf,
const double tn) {
359template <
int DIM,
typename BoundaryEleOp>
362 boost::shared_ptr<CommonData> common_data_ptr)
364 commonDataPtr(common_data_ptr) {}
366template <
int DIM,
typename BoundaryEleOp>
376 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
379 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
381 t_sum_t(
i) += alpha * t_traction(
i);
386 constexpr int ind[] = {0, 1, 2};
387 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
393template <
int DIM,
typename AssemblyBoundaryEleOp>
396 boost::shared_ptr<CommonData> common_data_ptr)
398 commonDataPtr(common_data_ptr) {}
400template <
int DIM,
typename AssemblyBoundaryEleOp>
403 EntitiesFieldData::EntData &data) {
411 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
413 auto &nf = AssemblyBoundaryEleOp::locF;
415 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
418 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
419 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
420 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
421 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
423 size_t nb_base_functions = data.getN().size2() / 3;
424 auto t_base = data.getFTensor1N<3>();
425 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
429 t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
431 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
432 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
435 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
437 auto ts_a = AssemblyBoundaryEleOp::getTStime();
439 auto sdf = surfaceDistanceFunction(
440 ts_a, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
441 t_total_traction(0), t_total_traction(1), t_total_traction(2));
443 auto t_grad_sdf = gradSurfaceDistanceFunction(
444 ts_a, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
445 t_total_traction(0), t_total_traction(1), t_total_traction(2));
447 auto tn = -t_traction(
i) * t_grad_sdf(
i);
451 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
453 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
462 t_cP(
i,
j) * t_disp(
j) +
463 c * (
sdf * t_grad_sdf(
i));
466 for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
467 const double beta = alpha * (t_base(
i) * t_normal(
i));
468 t_nf(
i) -= beta * t_rhs(
i);
473 for (; bb < nb_base_functions; ++bb)
486template <
int DIM,
typename AssemblyBoundaryEleOp>
489 const std::string col_field_name,
490 boost::shared_ptr<CommonData> common_data_ptr)
493 commonDataPtr(common_data_ptr) {
494 AssemblyBoundaryEleOp::sYmm =
false;
497template <
int DIM,
typename AssemblyBoundaryEleOp>
500 EntitiesFieldData::EntData &row_data,
501 EntitiesFieldData::EntData &col_data) {
508 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
509 auto &locMat = AssemblyBoundaryEleOp::locMat;
511 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
514 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
515 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
516 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
518 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
519 auto t_row_base = row_data.getFTensor1N<3>();
520 size_t nb_face_functions = row_data.getN().size2() / 3;
524 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
528 t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
530 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
533 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
535 auto ts_time = AssemblyBoundaryEleOp::getTStime();
537 auto sdf = surfaceDistanceFunction(
538 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
539 t_total_traction(0), t_total_traction(1), t_total_traction(2));
540 auto t_grad_sdf = gradSurfaceDistanceFunction(
541 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
542 t_total_traction(0), t_total_traction(1), t_total_traction(2));
544 auto tn = -t_traction(
i) * t_grad_sdf(
i);
548 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
550 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
553 t_res_dU(
i,
j) = kronecker_delta(
i,
j) + t_cP(
i,
j);
556 auto t_hess_sdf = hessSurfaceDistanceFunction(
557 ts_time, t_spatial_coords(0), t_spatial_coords(1),
558 t_spatial_coords(2), t_total_traction(0), t_total_traction(1),
559 t_total_traction(2));
562 (t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
563 t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
564 c *
sdf * t_hess_sdf(
i,
j);
568 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
570 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
572 const double row_base = t_row_base(
i) * t_normal(
i);
574 auto t_col_base = col_data.getFTensor0N(gg, 0);
575 for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
576 const double beta = alpha * row_base * t_col_base;
578 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
586 for (; rr < nb_face_functions; ++rr)
599template <
int DIM,
typename AssemblyBoundaryEleOp>
602 const std::string row_field_name,
const std::string col_field_name,
603 boost::shared_ptr<CommonData> common_data_ptr)
606 commonDataPtr(common_data_ptr) {
607 AssemblyBoundaryEleOp::sYmm =
false;
610template <
int DIM,
typename AssemblyBoundaryEleOp>
613 iNtegrate(EntitiesFieldData::EntData &row_data,
614 EntitiesFieldData::EntData &col_data) {
621 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
622 auto &locMat = AssemblyBoundaryEleOp::locMat;
624 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
627 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
628 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
629 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
631 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
632 auto t_row_base = row_data.getFTensor1N<3>();
633 size_t nb_face_functions = row_data.getN().size2() / 3;
635 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
639 t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
641 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
644 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
646 auto ts_time = AssemblyBoundaryEleOp::getTStime();
648 auto sdf = surfaceDistanceFunction(
649 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
650 t_total_traction(0), t_total_traction(1), t_total_traction(2));
651 auto t_grad_sdf = gradSurfaceDistanceFunction(
652 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
653 t_total_traction(0), t_total_traction(1), t_total_traction(2));
655 auto tn = -t_traction(
i) * t_grad_sdf(
i);
659 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
661 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
667 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
669 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
670 const double row_base = t_row_base(
i) * t_normal(
i);
672 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
673 for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
674 const double col_base = t_col_base(
i) * t_normal(
i);
675 const double beta = alpha * row_base * col_base;
677 t_mat(
i,
j) -= beta * t_res_dt(
i,
j);
685 for (; rr < nb_face_functions; ++rr)
698template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
700 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
701 std::string sigma, std::string u) {
703 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
705 using OpMixDivURhs =
typename B::template OpMixDivTimesU<3, DIM, DIM>;
706 using OpMixLambdaGradURhs =
typename B::template OpMixTensorTimesGradU<DIM>;
708 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
712 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
713 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
714 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
715 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
717 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
718 u, common_data_ptr->contactDispPtr()));
720 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
722 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
727 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(),
728 [](
double,
double,
double)
constexpr { return 1; }));
729 pip.push_back(
new OpMixLambdaGradURhs(sigma, mat_grad_ptr));
737 using OpMixLhs::OpMixLhs;
740 EntitiesFieldData::EntData &row_data,
741 EntitiesFieldData::EntData &col_data) {
743 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
744 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
746 if (side_fe_entity == side_fe_data) {
747 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
754template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEle>
757 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
758 std::string fe_domain_name, std::string sigma, std::string u,
759 std::string geom, ForcesAndSourcesCore::RuleHookFun rule) {
761 auto op_loop_side =
new OpLoopSide<DomainEle>(
762 m_field, fe_domain_name, DIM, Sev::noisy,
763 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
764 pip.push_back(op_loop_side);
766 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
769 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
772 using OpMixDivULhs =
typename B::template OpMixDivTimesVec<DIM>;
773 using OpLambdaGraULhs =
typename B::template OpMixTensorTimesGrad<DIM>;
777 auto unity = []() {
return 1; };
778 op_loop_side->getOpPtrVector().push_back(
779 new OpMixDivULhsSide(sigma, u, unity,
true));
780 op_loop_side->getOpPtrVector().push_back(
781 new OpLambdaGraULhsSide(sigma, u, unity,
true));
783 op_loop_side->getSideFEPtr()->getRuleHook = rule;
787template <
int DIM, AssemblyType A, IntegrationType I,
typename BoundaryEleOp>
789 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
790 std::string sigma, std::string u) {
795 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
797 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
798 u, common_data_ptr->contactDispPtr()));
799 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
800 sigma, common_data_ptr->contactTractionPtr()));
802 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
803 DIM, GAUSS>(sigma, u, common_data_ptr));
804 pip.push_back(
new typename C::template Assembly<A>::
805 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
806 sigma, sigma, common_data_ptr));
811template <
int DIM, AssemblyType A, IntegrationType I,
typename BoundaryEleOp>
813 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
814 std::string sigma, std::string u) {
819 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
821 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
822 u, common_data_ptr->contactDispPtr()));
823 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
824 sigma, common_data_ptr->contactTractionPtr()));
826 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
827 DIM, GAUSS>(sigma, common_data_ptr));
832template <
int DIM, IntegrationType I,
typename BoundaryEleOp>
834 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
840 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
841 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
842 sigma, common_data_ptr->contactTractionPtr()));
843 pip.push_back(
new typename C::template OpAssembleTotalContactTraction<DIM, I>(
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
def grad_sdf(t, x, y, z, tx, ty, tz)
def hess_sdf(t, x, y, z, tx, ty, tz)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr double t
plate stiffness
constexpr auto field_name
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
default operator for TRI element
double getMeasure()
get measure of element
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element