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

Inheritance diagram for ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >:

Collaboration diagram for ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >:

Public Member Functions
	OpAssembleTotalContactAreaImpl (boost::shared_ptr< CommonData > common_data_ptr, bool is_axisymmetric=false, boost::shared_ptr< Range > contact_range_ptr=nullptr)

MoFEMErrorCode	doWork (int side, EntityType type, EntData &data)

Public Attributes
SurfaceDistanceFunction	surfaceDistanceFunction = surface_distance_function

GradSurfaceDistanceFunction	gradSurfaceDistanceFunction

Private Attributes
boost::shared_ptr< CommonData >	commonDataPtr

bool	isAxisymmetric

boost::shared_ptr< Range >	contactRange

Detailed Description

template<int DIM, typename BoundaryEleOp>
struct ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >

Definition at line 460 of file ContactOps.hpp.

Constructor & Destructor Documentation

◆ OpAssembleTotalContactAreaImpl()

template<int DIM, typename BoundaryEleOp >

ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >::OpAssembleTotalContactAreaImpl	(	boost::shared_ptr< CommonData >	common_data_ptr,
		bool	is_axisymmetric = `false`,
		boost::shared_ptr< Range >	contact_range_ptr = `nullptr`
	)

Definition at line 653 of file ContactOps.hpp.

     : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
       commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric),
       contactRange(contact_range_ptr) {}

Member Function Documentation

◆ doWork()

template<int DIM, typename BoundaryEleOp >

MoFEMErrorCode ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >::doWork	(	int	side,
		EntityType	type,
		EntData &	data
	)

Examples: ContactOps.hpp.

Definition at line 662 of file ContactOps.hpp.

                                               {
   MoFEMFunctionBegin;
  
   auto fe_type = BoundaryEleOp::getFEType();
  
   const auto fe_ent = BoundaryEleOp::getFEEntityHandle();
  
   if (contactRange->find(fe_ent) != contactRange->end()) {
     FTensor::Index<'i', DIM> i;
     FTensor::Index<'j', DIM> j;
     FTensor::Tensor1<double, 2> t_sum_a{0., 0.};
  
     auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
     auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
     auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
  
     auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->contactDispGrad);
     auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
  
     const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
     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 = BoundaryEleOp::getTStime();
     auto ts_time_step = BoundaryEleOp::getTStimeStep();
     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 t_sdf = getFTensor0FromVec(v_sdf);
     auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
     for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
       double jacobian = 1.;
       if (isAxisymmetric) {
         jacobian = 2. * M_PI * t_coords(0); // Axisymmetric Jacobian
       }
       auto tn = -t_traction(i) * t_grad_sdf(i);
       auto c = constrain(t_sdf, tn);
       double alpha = t_w * jacobian;
  
       FTensor::Tensor2<double, DIM, DIM> F;
       FTensor::Tensor2<double, DIM, DIM> invF;
       FTensor::Tensor1<double, DIM> t_normal_current;
  
       F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
       auto det = determinantTensor(F);
       CHKERR invertTensor(F, det, invF);
       t_normal_current(i) = det * (invF(j, i) * t_normal_at_pts(j));
  
       alpha *= sqrt(t_normal_current(i) * t_normal_current(i));
  
       if (fe_type == MBTRI) {
         alpha /= 2;
       }
       if (c > 1e-12) {
         t_sum_a(0) += alpha; // real area
       }
       t_sum_a(1) += alpha; // Potential area
       ++t_w;
       ++t_traction;
       ++t_coords;
       ++t_sdf;
       ++t_grad_sdf;
  
       ++t_grad;
       ++t_normal_at_pts;
     }
     constexpr int ind[] = {3, 4};
     CHKERR VecSetValues(commonDataPtr->totalTraction, 2, ind, &t_sum_a(0),
                         ADD_VALUES);
   }
   MoFEMFunctionReturn(0);
 }