OpFaceSideMaterialForce Struct Reference

#include "users_modules/eshelbian_plasticity/src/EshelbianOperators.hpp"

Inheritance diagram for OpFaceSideMaterialForce:

Collaboration diagram for OpFaceSideMaterialForce:

Public Types
using	OP = VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator

Public Member Functions
	OpFaceSideMaterialForce (boost::shared_ptr< DataAtIntegrationPts > data_ptr)
MoFEMErrorCode	doWork (int side, EntityType type, EntData &data)
	Caluclate face material force and normal pressure at gauss points.

Private Attributes
boost::shared_ptr< DataAtIntegrationPts >	dataAtPts
	data at integration pts

Detailed Description

Definition at line 1004 of file EshelbianOperators.hpp.

Member Typedef Documentation

OP

using OpFaceSideMaterialForce::OP = VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator

Definition at line 1008 of file EshelbianOperators.hpp.

Constructor & Destructor Documentation

OpFaceSideMaterialForce()

OpFaceSideMaterialForce::OpFaceSideMaterialForce ( boost::shared_ptr< DataAtIntegrationPts >  data_ptr )

inline

Definition at line 1010 of file EshelbianOperators.hpp.

      : OP(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}

Member Function Documentation

doWork()

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

Caluclate face material force and normal pressure at gauss points.

Parameters

side
type
data

Returns

MoFEMErrorCode

Reconstruct the full gradient \(U=\nabla u\) on a surface from the symmetric part and the surface gradient.

Inputs:

• t_strain : \(\varepsilon=\tfrac12(U+U^\top)\) (symmetric strain on S),
• t_grad_u_gamma : \(u^\Gamma = U P\) (right-projected/surface gradient), with \(P=I-\mathbf N\otimes\mathbf N\),
• t_normal : (possibly non‑unit) surface normal.

Procedure (pointwise on S): 1) Normalize the normal \(\mathbf n=\mathbf N/\|\mathbf N\|\). 2) Form the residual \(R=\varepsilon-\operatorname{sym}(u^\Gamma)\), where \(\operatorname{sym}(A)=\tfrac12(A+A^\top)\). 3) Recover the normal directional derivative (a vector) \(\mathbf v=\partial_{\mathbf n}u=2R\mathbf n-(\mathbf n^\top R\,\mathbf n)\,\mathbf n\). 4) Assemble the full gradient \(U = u^\Gamma + \mathbf v\otimes \mathbf n\).

Properties (sanity checks):

• \(\tfrac12(U+U^\top)=\varepsilon\) (matches the given symmetric part),
• \(U P = u^\Gamma\) (tangential/right-projected columns unchanged),
• Only the normal column is updated via \(\mathbf v\otimes\mathbf n\).

Mapping to variables in this snippet:

• \(\varepsilon \leftrightarrow\) t_strain,
• \(u^\Gamma \leftrightarrow\) t_grad_u_gamma,
• \(\mathbf N \leftrightarrow\) t_normal (normalized into t_N),
• \(R \leftrightarrow\) t_R,
• \(U \leftrightarrow\) t_grad_u.

Precondition

t_normal is nonzero; t_strain is symmetric.

Note

All indices use Einstein summation; computation is local to the surface point.

Examples

/home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianOperators.cpp.

Definition at line 4413 of file EshelbianOperators.cpp.

                                                              {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
  FTENSOR_INDEX(SPACE_DIM, J);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, K);
  FTENSOR_INDEX(SPACE_DIM, M);
  FTENSOR_INDEX(SPACE_DIM, N);
  FTENSOR_INDEX(SPACE_DIM, L);
 
  const auto nb_gauss_pts = getGaussPts().size2();
  MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
      dataAtPts->faceMaterialForceAtPts, nb_gauss_pts);
  dataAtPts->normalPressureAtPts.resize(nb_gauss_pts, false);
  if (getNinTheLoop() == 0) {
    dataAtPts->faceMaterialForceAtPts.clear();
    dataAtPts->normalPressureAtPts.clear();
  }
  auto loop_size = getLoopSize();
  if (loop_size == 1) {
    auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
    auto pstatus = numebered_fe_ptr->getPStatus();
    if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
      loop_size = 2;
    }
  }
 
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
 
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_T = dataAtPts->getFTensorFaceMaterialForce(
      nb_gauss_pts); //< face material force
  auto t_p =
      getFTensor0FromVec(dataAtPts->normalPressureAtPts); //< normal pressure
  auto t_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
  auto t_u_gamma = dataAtPts->getFTensorSmallHybridDisp(nb_gauss_pts);
  auto t_grad_u_gamma = dataAtPts->getFTensorGradHybridDisp(nb_gauss_pts);
  auto t_strain = dataAtPts->getFTensorLogStretch(nb_gauss_pts);
  auto t_omega = dataAtPts->getFTensorRotAxis(nb_gauss_pts);
 
  FTensor::Tensor1<double, SPACE_DIM> t_N;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_A;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_R;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_grad_u;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_Proj;
 
  auto next = [&]() {
    ++t_normal;
    ++t_P;
    // ++t_grad_P;
    ++t_omega;
    ++t_u_gamma;
    ++t_grad_u_gamma;
    ++t_strain;
    ++t_T;
    ++t_p;
  };
 
  switch (EshelbianCore::energyReleaseSelector) {
  case GRIFFITH_FORCE:
    for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
      t_N(I) = t_normal(I);
      t_N.normalize();
 
      t_A(i, j) = levi_civita(i, j, k) * t_omega(k);
      t_R(i, k) = t_kd(i, k) + t_A(i, k);
      t_grad_u(i, j) = t_R(i, j) + t_strain(i, j);
 
      t_T(I) += t_N(J) * (t_grad_u(i, I) * t_P(i, J)) / loop_size;
      // note that works only for Hooke material, for nonlinear material we need
      // strain energy expressed by stress
      t_T(I) -= t_N(I) * ((t_strain(i, K) * t_P(i, K)) / 2.) / loop_size;
 
      t_p += t_N(I) *
             (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
             loop_size;
 
      next();
    }
    break;
  case GRIFFITH_SKELETON:
    for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
 
      // Normalize the normal
      t_N(I) = t_normal(I);
      t_N.normalize();
 
      // R = ε − sym(u^Γ)
      t_R(i, j) =
          t_strain(i, j) - 0.5 * (t_grad_u_gamma(i, j) + t_grad_u_gamma(j, i));
 
      // U = u^Γ + [2 R N − (Nᵀ R N) N] ⊗ N
      t_grad_u(i, J) =
          t_grad_u_gamma(i, J) +
          (2 * t_R(i, K) * t_N(K) - (t_R(k, L) * t_N(k) * t_N(L)) * t_N(i)) *
              t_N(J);
 
      t_T(I) += t_N(J) * (t_grad_u(i, I) * t_P(i, J)) / loop_size;
      // note that works only for Hooke material, for nonlinear material we need
      // strain energy expressed by stress
      t_T(I) -= t_N(I) * ((t_strain(i, K) * t_P(i, K)) / 2.) / loop_size;
 
      // calculate nominal face pressure
      t_p += t_N(I) *
             (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
             loop_size;
 
      next();
    }
    break;
 
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Grffith energy release "
            "selector not implemented");
  };
 
#ifndef NDEBUG
  auto side_fe_ptr = getSidePtrFE();
  auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
  auto pstatus = side_fe_mi_ptr->getPStatus();
  if (pstatus) {
    auto owner = side_fe_mi_ptr->getOwnerProc();
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
        << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
        << getNinTheLoop() << " loop size " << getLoopSize();
  }
#endif // NDEBUG
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

dataAtPts

boost::shared_ptr<DataAtIntegrationPts> OpFaceSideMaterialForce::dataAtPts

private

data at integration pts

Definition at line 1017 of file EshelbianOperators.hpp.

---

The documentation for this struct was generated from the following files:

• users_modules/eshelbian_plasticity/src/EshelbianOperators.hpp
• users_modules/eshelbian_plasticity/src/impl/EshelbianOperators.cpp