#include <tutorials/adv-1/src/ContactOps.hpp>

Inheritance diagram for ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >:

Collaboration diagram for ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >:

Public Member Functions
	OpConstrainBoundaryLhs_dUImpl (const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, bool is_axisymmetric=false)

MoFEMErrorCode	iNtegrate (EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)

Public Attributes
SurfaceDistanceFunction	surfaceDistanceFunction = surface_distance_function

GradSurfaceDistanceFunction	gradSurfaceDistanceFunction

HessSurfaceDistanceFunction	hessSurfaceDistanceFunction

boost::shared_ptr< CommonData >	commonDataPtr

bool	isAxisymmetric

Detailed Description

template<int DIM, typename AssemblyBoundaryEleOp>
struct ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >

Definition at line 511 of file ContactOps.hpp.

Constructor & Destructor Documentation

◆ OpConstrainBoundaryLhs_dUImpl()

template<int DIM, typename AssemblyBoundaryEleOp >

ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::OpConstrainBoundaryLhs_dUImpl	(	const std::string	row_field_name,
		const std::string	col_field_name,
		boost::shared_ptr< CommonData >	common_data_ptr,
		bool	is_axisymmetric = `false`
	)

Definition at line 950 of file ContactOps.hpp.

     : AssemblyBoundaryEleOp(row_field_name, col_field_name,
                             AssemblyBoundaryEleOp::OPROWCOL),
       commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
   AssemblyBoundaryEleOp::sYmm = false;
 }

Member Function Documentation

◆ iNtegrate()

template<int DIM, typename AssemblyBoundaryEleOp >

MoFEMErrorCode ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::iNtegrate	(	EntitiesFieldData::EntData &	row_data,
		EntitiesFieldData::EntData &	col_data
	)

Examples: ContactOps.hpp.

Definition at line 962 of file ContactOps.hpp.

                                         {
   MoFEMFunctionBegin;
  
   FTensor::Index<'i', DIM> i;
   FTensor::Index<'j', DIM> j;
   FTensor::Index<'k', DIM> k;
  
   const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
   auto &locMat = AssemblyBoundaryEleOp::locMat;
  
   auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
   auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
   auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
  
   auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
   auto t_row_base = row_data.getFTensor1N<3>();
   size_t nb_face_functions = row_data.getN().size2() / 3;
  
   constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  
   auto m_spatial_coords = get_spatial_coords(
       BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
       getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
   auto m_normals_at_pts = get_normalize_normals(
       BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
  
   auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
  
   auto ts_time = AssemblyBoundaryEleOp::getTStime();
   auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
  
   // placeholder to pass boundary block id to python
   int block_id = 0;
  
   auto v_sdf =
       surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                               m_spatial_coords, m_normals_at_pts, block_id);
  
   auto m_grad_sdf =
       gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                   m_spatial_coords, m_normals_at_pts, block_id);
  
   auto m_hess_sdf =
       hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                   m_spatial_coords, m_normals_at_pts, block_id);
  
   auto t_sdf = getFTensor0FromVec(v_sdf);
   auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
   auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
  
   for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
  
     double jacobian = 1.;
     if (isAxisymmetric) {
       jacobian = 2. * M_PI * t_coords(0);
     }
     const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
  
     auto tn = -t_traction(i) * t_grad_sdf(i);
     auto c = constrain(t_sdf, tn);
  
     FTensor::Tensor2<double, DIM, DIM> t_cP;
     t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
     FTensor::Tensor2<double, DIM, DIM> t_cQ;
     t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
  
     FTensor::Tensor2<double, DIM, DIM> t_res_dU;
     t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
  
     if (c > 0) {
       t_res_dU(i, j) +=
           (c * cn_contact) *
               (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
                t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
           c * t_sdf * t_hess_sdf(i, j);
     }
  
     size_t rr = 0;
     for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
  
       auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
  
       const double row_base = t_row_base(i) * t_normal(i);
  
       auto t_col_base = col_data.getFTensor0N(gg, 0);
       for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
         const double beta = alpha * row_base * t_col_base;
  
         t_mat(i, j) -= beta * t_res_dU(i, j);
  
         ++t_col_base;
         ++t_mat;
       }
  
       ++t_row_base;
     }
     for (; rr < nb_face_functions; ++rr)
       ++t_row_base;
  
     ++t_traction;
     ++t_coords;
     ++t_w;
     ++t_normal;
     ++t_sdf;
     ++t_grad_sdf;
     ++t_hess_sdf;
   }
  
   MoFEMFunctionReturn(0);
 }