EshelbianPlasticity Namespace Reference

Classes
struct	AddHOOps
struct	AddHOOps< 2, 3, 3 >
struct	AddHOOps< 3, 3, 3 >
struct	BcDisp
struct	BcRot
struct	CGGUserPolynomialBase
struct	ContactTree
struct	DataAtIntegrationPts
struct	EshelbianCore
struct	FaceRule
struct	FaceUserDataOperatorStabAssembly
struct	HMHHencky
struct	HMHNeohookean
struct	HMHPMooneyRivlinWriggersEq63
struct	HMHStVenantKirchhoff
struct	isEq
struct	OpAssembleBasic
struct	OpAssembleFace
struct	OpAssembleVolume
struct	OpBrokenTractionBc
struct	OpCalculateEshelbyStress
struct	OpCalculateRotationAndSpatialGradient
struct	OpConstrainBoundaryHDivLhs_dU
struct	OpConstrainBoundaryHDivLhs_dU< A, IntegrationType::GAUSS >
struct	OpConstrainBoundaryHDivRhs
struct	OpConstrainBoundaryHDivRhs< A, IntegrationType::GAUSS >
struct	OpConstrainBoundaryL2Lhs_dP
struct	OpConstrainBoundaryL2Lhs_dP< A, IntegrationType::GAUSS >
struct	OpConstrainBoundaryL2Lhs_dU
struct	OpConstrainBoundaryL2Rhs
struct	OpDispBc
struct	OpDispBcImpl
struct	OpDispBcImpl< A, GAUSS >
struct	OpHMHH
struct	OpJacobian
struct	OpMoveNode
struct	OpPostProcDataStructure
struct	OpRotationBc
struct	OpRotationBcImpl
struct	OpRotationBcImpl< A, GAUSS >
struct	OpSpatialConsistency_dBubble_dBubble
struct	OpSpatialConsistency_dBubble_domega
struct	OpSpatialConsistency_dBubble_dP
struct	OpSpatialConsistency_dP_domega
struct	OpSpatialConsistency_dP_dP
struct	OpSpatialConsistencyBubble
struct	OpSpatialConsistencyDivTerm
struct	OpSpatialConsistencyP
struct	OpSpatialEquilibrium
struct	OpSpatialEquilibrium_dw_dP
struct	OpSpatialEquilibrium_dw_dw
struct	OpSpatialPhysical
struct	OpSpatialPhysical_du_dBubble
struct	OpSpatialPhysical_du_domega
struct	OpSpatialPhysical_du_dP
struct	OpSpatialPhysical_du_du
struct	OpSpatialPrj
struct	OpSpatialPrj_dx_dw
struct	OpSpatialPrj_dx_dx
struct	OpSpatialRotation
struct	OpSpatialRotation_domega_dBubble
struct	OpSpatialRotation_domega_domega
struct	OpSpatialRotation_domega_dP
struct	OpTreeSearch
struct	PhysicalEquations
struct	SetIntegrationAtFrontVolume
struct	TensorTypeExtractor
struct	TensorTypeExtractor< FTensor::PackPtr< T, I > >
struct	TractionBc
struct	VolRule
	Set integration rule on element. More...
struct	VolUserDataOperatorStabAssembly

Typedefs
using	MatrixPtr = boost::shared_ptr< MatrixDouble >
using	VectorPtr = boost::shared_ptr< VectorDouble >
using	EntData = EntitiesFieldData::EntData
using	UserDataOperator = ForcesAndSourcesCore::UserDataOperator
using	VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
using	FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator
typedef std::vector< BcDisp >	BcDispVec
typedef std::vector< BcRot >	BcRotVec
typedef std::vector< Range >	TractionFreeBc
typedef std::vector< TractionBc >	TractionBcVec

Enumerations
enum	RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION }
enum	StretchSelector { LINEAR, LOG }
enum	MultiPointRhsType { U, P }

Functions
MoFEMErrorCode	CGG_BubbleBase_MBTET (const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
	Calculate CGGT tonsorial bubble base. More...
auto	diff_deviator (FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
auto	diff_tensor ()
auto	symm_L_tensor ()
auto	diff_symmetrize ()
auto	checkSdf (EntityHandle fe_ent, std::map< int, Range > &sdf_map_range)
template<typename OP_PTR >
auto	getSdf (OP_PTR op_ptr, MatrixDouble &contact_disp, int block_id, bool eval_hessian)
template<typename T1 >
auto	multiPoint (std::array< double, 3 > &unit_ray, std::array< double, 3 > &point, std::array< double, 9 > &elem_point_nodes, std::array< double, 9 > &elem_traction_nodes, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
	Calculate points data on contact surfaces. More...
template<typename T1 >
auto	multiMasterPoint (ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
template<typename T1 >
auto	multiSlavePoint (ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
template<typename T1 >
auto	multiGetGap (ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
template<typename T1 , typename T2 , typename T3 >
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)
template<typename T >
auto	get_rotation_form_vector (FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
template<typename T >
auto	get_diff_rotation_form_vector (FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
template<typename T >
auto	get_diff2_rotation_form_vector (FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
auto	is_eq (const double &a, const double &b)
template<int DIM>
auto	get_uniq_nb (double *ptr)
template<int DIM>
auto	sort_eigen_vals (FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)

Variables
constexpr static auto	size_symm = (3 * (3 + 1)) / 2

Typedef Documentation

BcDispVec

typedef std::vector<BcDisp> EshelbianPlasticity::BcDispVec

Definition at line 358 of file EshelbianPlasticity.hpp.

BcRotVec

typedef std::vector<BcRot> EshelbianPlasticity::BcRotVec

Definition at line 367 of file EshelbianPlasticity.hpp.

EntData

using EshelbianPlasticity::EntData = typedef EntitiesFieldData::EntData

Definition at line 49 of file EshelbianPlasticity.hpp.

FaceUserDataOperator

using EshelbianPlasticity::FaceUserDataOperator = typedef FaceElementForcesAndSourcesCore::UserDataOperator

Examples

NavierStokesElement.hpp.

Definition at line 52 of file EshelbianPlasticity.hpp.

MatrixPtr

using EshelbianPlasticity::MatrixPtr = typedef boost::shared_ptr<MatrixDouble>

Definition at line 46 of file EshelbianPlasticity.hpp.

TractionBcVec

typedef std::vector<TractionBc> EshelbianPlasticity::TractionBcVec

Definition at line 378 of file EshelbianPlasticity.hpp.

TractionFreeBc

typedef std::vector<Range> EshelbianPlasticity::TractionFreeBc

Definition at line 369 of file EshelbianPlasticity.hpp.

UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

Definition at line 50 of file EshelbianPlasticity.hpp.

VectorPtr

using EshelbianPlasticity::VectorPtr = typedef boost::shared_ptr<VectorDouble>

Definition at line 47 of file EshelbianPlasticity.hpp.

VolUserDataOperator

using EshelbianPlasticity::VolUserDataOperator = typedef VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 51 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

MultiPointRhsType

enum EshelbianPlasticity::MultiPointRhsType

Enumerator
U
P

Definition at line 197 of file EshelbianContact.cpp.

{ U, P };

RotSelector

enum EshelbianPlasticity::RotSelector

Enumerator
SMALL_ROT
MODERATE_ROT
LARGE_ROT
NO_H1_CONFIGURATION

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 43 of file EshelbianPlasticity.hpp.

{ SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION };

StretchSelector

enum EshelbianPlasticity::StretchSelector

Enumerator
LINEAR
LOG

Examples

EshelbianPlasticity.cpp.

Definition at line 44 of file EshelbianPlasticity.hpp.

{ LINEAR, LOG };

Function Documentation

CGG_BubbleBase_MBTET()

MoFEMErrorCode EshelbianPlasticity::CGG_BubbleBase_MBTET	(	const int	p,
		const double *	N,
		const double *	diffN,
		FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &	phi,
		const int	gdim
	)

Calculate CGGT tonsorial bubble base.

See details in [20]

Parameters

p	polynomial order
N	shape functions
diffN	direvatives of shape functions
l2_base	base functions for l2 space
diff_l2_base	direvatives base functions for l2 space
phi	returned base functions
gdim	number of integration points

Returns

MoFEMErrorCode

Examples

EshelbianPlasticity.cpp.

Definition at line 20 of file CGGTonsorialBubbleBase.cpp.

                                                    {
  MoFEMFunctionBeginHot;
 
  // FIXME: In this implementation symmetry of mix derivatives is not exploited
 
  if (p < 1)
    MoFEMFunctionReturnHot(0);
 
  Index<'i', 3> i;
  Index<'j', 3> j;
  Index<'k', 3> k;
  Index<'m', 3> m;
  Index<'f', 3> f;
 
  Tensor1<double, 3> t_diff_n[4];
  {
    Tensor1<PackPtr<const double *, 3>, 3> t_diff_n_tmp(&diffN[0], &diffN[1],
                                                        &diffN[2]);
    for (int ii = 0; ii != 4; ++ii) {
      t_diff_n[ii](i) = t_diff_n_tmp(i);
      ++t_diff_n_tmp;
    }
  }
 
  Tensor1<double, 3> t_diff_ksi[3];
  for (int ii = 0; ii != 3; ++ii)
    t_diff_ksi[ii](i) = t_diff_n[ii + 1](i) - t_diff_n[0](i);
 
  int lp = p >= 2 ? p - 2 + 1 : 0;
  VectorDouble l[3] = {VectorDouble(lp + 1), VectorDouble(lp + 1),
                       VectorDouble(lp + 1)};
  MatrixDouble diff_l[3] = {MatrixDouble(3, lp + 1), MatrixDouble(3, lp + 1),
                            MatrixDouble(3, lp + 1)};
  MatrixDouble diff2_l[3] = {MatrixDouble(9, lp + 1), MatrixDouble(9, lp + 1),
                             MatrixDouble(9, lp + 1)};
 
  for (int ii = 0; ii != 3; ++ii)
    diff2_l[ii].clear();
 
  for (int gg = 0; gg != gdim; ++gg) {
 
    const int node_shift = gg * 4;
 
    for (int ii = 0; ii != 3; ++ii) {
 
      auto &t_diff_ksi_ii = t_diff_ksi[ii];
      auto &l_ii = l[ii];
      auto &diff_l_ii = diff_l[ii];
      auto &diff2_l_ii = diff2_l[ii];
 
      double ksi_ii = N[node_shift + ii + 1] - N[node_shift + 0];
 
      CHKERR Legendre_polynomials(lp, ksi_ii, &t_diff_ksi_ii(0),
                                  &*l_ii.data().begin(),
                                  &*diff_l_ii.data().begin(), 3);
 
      for (int l = 1; l < lp; ++l) {
        const double a = ((2 * (double)l + 1) / ((double)l + 1));
        const double b = ((double)l / ((double)l + 1));
        for (int d0 = 0; d0 != 3; ++d0)
          for (int d1 = 0; d1 != 3; ++d1) {
            const int r = 3 * d0 + d1;
            diff2_l_ii(r, l + 1) = a * (t_diff_ksi_ii(d0) * diff_l_ii(d1, l) +
                                        t_diff_ksi_ii(d1) * diff_l_ii(d0, l) +
                                        ksi_ii * diff2_l_ii(r, l)) -
                                   b * diff2_l_ii(r, l - 1);
          }
      }
    }
 
    const double n[] = {N[node_shift + 0], N[node_shift + 1], N[node_shift + 2],
                        N[node_shift + 3]};
 
    Tensor2<double, 3, 3> t_bk;
    Tensor3<double, 3, 3, 3> t_bk_diff;
    const int tab[4][4] = {
        {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
    t_bk(i, j) = 0;
    t_bk_diff(i, j, k) = 0;
    for (int ii = 0; ii != 3; ++ii) {
      const int i0 = tab[ii][0];
      const int i1 = tab[ii][1];
      const int i2 = tab[ii][2];
      const int i3 = tab[ii][3];
      auto &t_diff_n_i0 = t_diff_n[i0];
      auto &t_diff_n_i1 = t_diff_n[i1];
      auto &t_diff_n_i2 = t_diff_n[i2];
      auto &t_diff_n_i3 = t_diff_n[i3];
      Tensor2<double, 3, 3> t_k;
      t_k(i, j) = t_diff_n_i3(i) * t_diff_n_i3(j);
      const double b = n[i0] * n[i1] * n[i2];
      t_bk(i, j) += b * t_k(i, j);
      Tensor1<double, 3> t_diff_b;
      t_diff_b(i) = t_diff_n_i0(i) * n[i1] * n[i2] +
                    t_diff_n_i1(i) * n[i0] * n[i2] +
                    t_diff_n_i2(i) * n[i0] * n[i1];
      t_bk_diff(i, j, k) += t_k(i, j) * t_diff_b(k);
    }
 
    int zz = 0;
    for (int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
 
      for (int ii = 0; ii <= o; ++ii)
        for (int jj = 0; (ii + jj) <= o; ++jj) {
 
          const int kk = o - ii - jj;
 
          auto get_diff_l = [&](const int y, const int i) {
            return Tensor1<double, 3>(diff_l[y](0, i), diff_l[y](1, i),
                                      diff_l[y](2, i));
          };
          auto get_diff2_l = [&](const int y, const int i) {
            return Tensor2<double, 3, 3>(
                diff2_l[y](0, i), diff2_l[y](1, i), diff2_l[y](2, i),
                diff2_l[y](3, i), diff2_l[y](4, i), diff2_l[y](5, i),
                diff2_l[y](6, i), diff2_l[y](7, i), diff2_l[y](8, i));
          };
 
          auto l_i = l[0][ii];
          auto t_diff_i = get_diff_l(0, ii);
          auto t_diff2_i = get_diff2_l(0, ii);
          auto l_j = l[1][jj];
          auto t_diff_j = get_diff_l(1, jj);
          auto t_diff2_j = get_diff2_l(1, jj);
          auto l_k = l[2][kk];
          auto t_diff_k = get_diff_l(2, kk);
          auto t_diff2_k = get_diff2_l(2, kk);
 
          Tensor1<double, 3> t_diff_l2;
          t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
                         t_diff_k(i) * l_i * l_j;
          Tensor2<double, 3, 3> t_diff2_l2;
          t_diff2_l2(i, j) =
              t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
              t_diff_i(i) * l_j * t_diff_k(j) +
 
              t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
              t_diff_j(i) * l_i * t_diff_k(j) +
 
              t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
              t_diff_k(i) * l_i * t_diff_j(j);
 
          for (int dd = 0; dd != 3; ++dd) {
 
            Tensor2<double, 3, 3> t_axial_diff;
            t_axial_diff(i, j) = 0;
            for (int mm = 0; mm != 3; ++mm)
              t_axial_diff(dd, mm) = t_diff_l2(mm);
 
            Tensor3<double, 3, 3, 3> t_A_diff;
            t_A_diff(i, j, k) = levi_civita(i, j, m) * t_axial_diff(m, k);
            Tensor2<double, 3, 3> t_curl_A;
            t_curl_A(i, j) = levi_civita(j, m, f) * t_A_diff(i, f, m);
            Tensor3<double, 3, 3, 3> t_curl_A_bK_diff;
            t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
 
            Tensor3<double, 3, 3, 3> t_axial_diff2;
            t_axial_diff2(i, j, k) = 0;
            for (int mm = 0; mm != 3; ++mm)
              for (int nn = 0; nn != 3; ++nn)
                t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
            Tensor4<double, 3, 3, 3, 3> t_A_diff2;
            t_A_diff2(i, j, k, f) =
                levi_civita(i, j, m) * t_axial_diff2(m, k, f);
            Tensor3<double, 3, 3, 3> t_curl_A_diff2;
            t_curl_A_diff2(i, j, k) =
                levi_civita(j, m, f) * t_A_diff2(i, f, m, k);
            Tensor3<double, 3, 3, 3> t_curl_A_diff2_bK;
            t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
 
            t_phi(i, j) = levi_civita(j, m, f) * (t_curl_A_bK_diff(i, f, m) +
                                                  t_curl_A_diff2_bK(i, f, m));
 
            ++t_phi;
            ++zz;
          }
        }
    }
    if (zz != NBVOLUMETET_CCG_BUBBLE(p))
      SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "Wrong number of base functions %d != %d", zz,
               NBVOLUMETET_CCG_BUBBLE(p));
  }
 
  MoFEMFunctionReturnHot(0);
}

checkSdf()

auto EshelbianPlasticity::checkSdf	(	EntityHandle	fe_ent,
		std::map< int, Range > &	sdf_map_range
	)

Definition at line 17 of file EshelbianContact.cpp.

                                                                      {
  for (auto &m_sdf : sdf_map_range) {
    if (m_sdf.second.find(fe_ent) != m_sdf.second.end())
      return m_sdf.first;
  }
  return -1;
}

diff_deviator()

auto EshelbianPlasticity::diff_deviator ( FTensor::Ddg< double, 3, 3 > &&  t_diff_stress )

inline

Definition at line 12 of file EshelbianAux.hpp.

                                                                    {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  FTensor::Ddg<double, 3, 3> t_diff_deviator;
  t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
 
  constexpr double third = boost::math::constants::third<double>();
 
  t_diff_deviator(0, 0, 0, 0) -= third;
  t_diff_deviator(0, 0, 1, 1) -= third;
 
  t_diff_deviator(1, 1, 0, 0) -= third;
  t_diff_deviator(1, 1, 1, 1) -= third;
 
  t_diff_deviator(2, 2, 0, 0) -= third;
  t_diff_deviator(2, 2, 1, 1) -= third;
 
  t_diff_deviator(0, 0, 2, 2) -= third;
  t_diff_deviator(1, 1, 2, 2) -= third;
  t_diff_deviator(2, 2, 2, 2) -= third;
 
  return t_diff_deviator;
}

diff_symmetrize()

auto EshelbianPlasticity::diff_symmetrize ( )

inline

Definition at line 67 of file EshelbianAux.hpp.

                              {
 
  FTensor::Index<'i', SPACE_DIM> i;
  FTensor::Index<'j', SPACE_DIM> j;
  FTensor::Index<'k', SPACE_DIM> k;
  FTensor::Index<'l', SPACE_DIM> l;
 
  FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_diff;
 
  t_diff(i, j, k, l) = 0;
  t_diff(0, 0, 0, 0) = 1;
  t_diff(1, 1, 1, 1) = 1;
 
  t_diff(1, 0, 1, 0) = 0.5;
  t_diff(1, 0, 0, 1) = 0.5;
 
  t_diff(0, 1, 0, 1) = 0.5;
  t_diff(0, 1, 1, 0) = 0.5;
 
  if constexpr (SPACE_DIM == 3) {
    t_diff(2, 2, 2, 2) = 1;
 
    t_diff(2, 0, 2, 0) = 0.5;
    t_diff(2, 0, 0, 2) = 0.5;
    t_diff(0, 2, 0, 2) = 0.5;
    t_diff(0, 2, 2, 0) = 0.5;
 
    t_diff(2, 1, 2, 1) = 0.5;
    t_diff(2, 1, 1, 2) = 0.5;
    t_diff(1, 2, 1, 2) = 0.5;
    t_diff(1, 2, 2, 1) = 0.5;
  }
 
  return t_diff;
}

diff_tensor()

auto EshelbianPlasticity::diff_tensor ( )

inline

Definition at line 41 of file EshelbianAux.hpp.

                          {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Ddg<double, 3, 3> t_diff;
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
  return t_diff;
};

get_diff2_rotation_form_vector()

template<typename T >

auto EshelbianPlasticity::get_diff2_rotation_form_vector	(	FTensor::Tensor1< T, 3 > &	t_omega,
		RotSelector	rotSelector = LARGE_ROT
	)

Examples

EshelbianOperators.cpp.

Definition at line 108 of file EshelbianOperators.cpp.

                                                                         {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Tensor4<D, 3, 3, 3, 3> t_diff2_R;
  FTensor::Index<'i', 3> i;
 
  constexpr double eps = 1e-10;
  for (int l = 0; l != 3; ++l) {
    FTensor::Tensor1<std::complex<double>, 3> t_omega_c;
    t_omega_c(i) = t_omega(i);
    t_omega_c(l) += std::complex<double>(0, eps);
    auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
    for (int i = 0; i != 3; ++i) {
      for (int j = 0; j != 3; ++j) {
        for (int k = 0; k != 3; ++k) {
          t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
        }
      }
    }
  }
 
  return t_diff2_R;
}

get_diff_rotation_form_vector()

template<typename T >

auto EshelbianPlasticity::get_diff_rotation_form_vector	(	FTensor::Tensor1< T, 3 > &	t_omega,
		RotSelector	rotSelector = LARGE_ROT
	)

Examples

EshelbianOperators.cpp.

Definition at line 60 of file EshelbianOperators.cpp.

                                                                        {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  FTensor::Tensor3<D, 3, 3, 3> t_diff_R;
 
  const auto angle =
      sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
  if (rotSelector == SMALL_ROT) {
    t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
    return t_diff_R;
  }
 
  const auto ss = sin(angle);
  const auto a = ss / angle;
  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
 
  FTensor::Tensor2<D, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
 
  const auto angle2 = angle * angle;
  const auto cc = cos(angle);
  const auto diff_a = (angle * cc - ss) / angle2;
  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
  if (rotSelector == MODERATE_ROT)
    return t_diff_R;
 
  const auto ss_2 = sin(angle / 2.);
  const auto cc_2 = cos(angle / 2.);
  const auto b = 2. * ss_2 * ss_2 / angle2;
  t_diff_R(i, j, k) +=
      b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
           FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
  const auto diff_b =
      (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
  t_diff_R(i, j, k) +=
      diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
 
  return t_diff_R;
}

get_rotation_form_vector()

template<typename T >

auto EshelbianPlasticity::get_rotation_form_vector	(	FTensor::Tensor1< T, 3 > &	t_omega,
		RotSelector	rotSelector = LARGE_ROT
	)

Examples

EshelbianOperators.cpp.

Definition at line 24 of file EshelbianOperators.cpp.

                                                                   {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Tensor2<D, 3, 3> t_R;
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  t_R(i, j) = t_kd(i, j);
 
  FTensor::Tensor2<D, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
 
  if (rotSelector == SMALL_ROT) {
    t_R(i, j) += t_Omega(i, j);
    return t_R;
  }
 
  const auto angle =
      sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
  const auto a = sin(angle) / angle;
  t_R(i, j) += a * t_Omega(i, j);
  if (rotSelector == MODERATE_ROT)
    return t_R;
 
  const auto ss_2 = sin(angle / 2.) / angle;
  const auto b = 2. * ss_2 * ss_2;
  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
 
  return t_R;
}

get_uniq_nb()

template<int DIM>

auto EshelbianPlasticity::get_uniq_nb ( double *  ptr )

inline

Examples

EshelbianOperators.cpp.

Definition at line 148 of file EshelbianOperators.cpp.

                                                        {
  std::array<double, DIM> tmp;
  std::copy(ptr, &ptr[DIM], tmp.begin());
  std::sort(tmp.begin(), tmp.end());
  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
  constexpr double eps = std::numeric_limits<float>::epsilon();
  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
}

getSdf()

template<typename OP_PTR >

auto EshelbianPlasticity::getSdf	(	OP_PTR	op_ptr,
		MatrixDouble &	contact_disp,
		int	block_id,
		bool	eval_hessian
	)

Definition at line 26 of file EshelbianContact.cpp.

                               {
 
  auto nb_gauss_pts = op_ptr->getGaussPts().size2();
 
  auto ts_time = op_ptr->getTStime();
  auto ts_time_step = op_ptr->getTStimeStep();
 
  auto m_spatial_coords = ContactOps::get_spatial_coords(
      op_ptr->getFTensor1CoordsAtGaussPts(),
      getFTensor1FromMat<3>(contact_disp), nb_gauss_pts);
  auto m_normals_at_pts = ContactOps::get_normalize_normals(
      op_ptr->getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
  auto v_sdf = ContactOps::surface_distance_function(
      ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
      block_id);
  auto m_grad_sdf = ContactOps::grad_surface_distance_function(
      ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
      block_id);
 
#ifndef NDEBUG
  if (v_sdf.size() != nb_gauss_pts)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong number of integration pts");
  if (m_grad_sdf.size2() != nb_gauss_pts)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong number of integration pts");
  if (m_grad_sdf.size1() != 3)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Should be size of 3");
#endif // NDEBUG
 
  if (eval_hessian) {
    auto m_hess_sdf = ContactOps::hess_surface_distance_function(
        ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
        block_id);
    return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
                           m_hess_sdf);
  } else {
 
    return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
                           MatrixDouble(6, nb_gauss_pts, 0.));
  }
};

is_eq()

auto EshelbianPlasticity::is_eq ( const double &  a, const double &  b  )

inline

Examples

EshelbianOperators.cpp.

Definition at line 144 of file EshelbianOperators.cpp.

                                                    {
  return isEq::check(a, b);
};

multiGetGap()

template<typename T1 >

auto EshelbianPlasticity::multiGetGap	(	ContactTree::FaceData *	face_data_ptr,
		FTensor::Tensor1< T1, 3 > &	t_spatial_coords
	)

Evaluate gap and tractions between master and slave points

Definition at line 173 of file EshelbianContact.cpp.

                                                          {
 
  auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
      multiSlavePoint(face_data_ptr, t_spatial_coords);
  auto [t_master_point_current, t_master_traction_current, t_master_normal] =
      multiMasterPoint(face_data_ptr, t_spatial_coords);
 
  FTensor::Index<'i', 3> i;
 
  FTensor::Tensor1<double, 3> t_normal;
  t_normal(i) = t_master_normal(i) - t_slave_normal(i);
  t_normal.normalize();
 
  auto gap = t_normal(i) * (t_slave_point_current(i) - t_spatial_coords(i));
  auto tn_master = t_master_traction_current(i) * t_normal(i);
  auto tn_slave = t_slave_traction_current(i) * t_normal(i);
  auto tn = std::max(-tn_master, tn_slave);
 
  return std::make_tuple(gap, tn_master, tn_slave,
                         ContactOps::constrain(gap, tn),
                         t_master_traction_current, t_slave_traction_current);
};

multiMasterPoint()

template<typename T1 >

auto EshelbianPlasticity::multiMasterPoint	(	ContactTree::FaceData *	face_data_ptr,
		FTensor::Tensor1< T1, 3 > &	t_spatial_coords
	)

Definition at line 144 of file EshelbianContact.cpp.

  {
 
  return multiPoint(face_data_ptr->unitRay, face_data_ptr->masterPoint,
                    face_data_ptr->masterPointNodes,
                    face_data_ptr->masterTractionNodes, t_spatial_coords);
};

multiPoint()

template<typename T1 >

auto EshelbianPlasticity::multiPoint	(	std::array< double, 3 > &	unit_ray,
		std::array< double, 3 > &	point,
		std::array< double, 9 > &	elem_point_nodes,
		std::array< double, 9 > &	elem_traction_nodes,
		FTensor::Tensor1< T1, 3 > &	t_spatial_coords
	)

Calculate points data on contact surfaces.

Template Parameters

T1

Parameters

unit_ray
point
elem_point_nodes
elem_traction_nodes
t_spatial_coords

Returns

auto

Definition at line 82 of file EshelbianContact.cpp.

                                             {
  FTensor::Index<'i', 3> i;
 
  auto t_unit_ray = getFTensor1FromPtr<3>(unit_ray.data());
  auto t_point = getFTensor1FromPtr<3>(point.data());
 
  auto get_normal = [](auto &ele_coords) {
    FTensor::Tensor1<double, 3> t_normal;
    Tools::getTriNormal(ele_coords.data(), &t_normal(0));
    return t_normal;
  };
 
  auto t_normal = get_normal(elem_point_nodes);
  t_normal(i) /= t_normal.l2();
 
  auto sn = t_normal(i) * t_point(i);
  auto nm = t_normal(i) * t_spatial_coords(i);
  auto nr = t_normal(i) * t_unit_ray(i);
 
  auto gamma = (sn - nm) / nr;
 
  FTensor::Tensor1<T1, 3> t_point_current;
  t_point_current(i) = t_spatial_coords(i) + gamma * t_unit_ray(i);
 
  auto get_local_point_shape_functions = [&](auto &&t_elem_coords,
                                             auto &t_point) {
    std::array<T1, 2> loc_coords;
    CHK_THROW_MESSAGE(
        Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
            &t_elem_coords(0, 0), &t_point(0), 1, loc_coords.data()),
        "get local coords");
    return FTensor::Tensor1<T1, 3>{N_MBTRI0(loc_coords[0], loc_coords[1]),
                                   N_MBTRI1(loc_coords[0], loc_coords[1]),
                                   N_MBTRI2(loc_coords[0], loc_coords[1])};
  };
 
  auto eval_position = [&](auto &&t_field, auto &t_shape_fun) {
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Tensor1<T1, 3> t_point_out;
    t_point_out(i) = t_shape_fun(j) * t_field(j, i);
    return t_point_out;
  };
 
  auto t_shape_fun = get_local_point_shape_functions(
      getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_point_current);
  auto t_slave_point_updated = eval_position(
      getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_shape_fun);
  auto t_traction_updated = eval_position(
      getFTensor2FromPtr<3, 3>(elem_traction_nodes.data()), t_shape_fun);
 
  return std::make_tuple(t_slave_point_updated, t_traction_updated, t_normal);
};

multiPointRhs()

template<typename T1 , typename T2 , typename T3 >

auto EshelbianPlasticity::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
	)

Caluclate rhs term for contact

Definition at line 202 of file EshelbianContact.cpp.

  {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  FTensor::Tensor1<T2, 3> t_u;
  t_u(i) = t_spatial_coords(i) - t_coords(i);
 
  FTensor::Tensor1<typename FTensor::promote<T2, T3>::V, 3> t_rhs;
 
  switch (type) {
  case MultiPointRhsType::U: {
 
    auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
        multiSlavePoint(face_data_ptr, t_spatial_coords);
    auto [t_master_point_current, t_master_traction_current, t_master_normal] =
        multiMasterPoint(face_data_ptr, t_spatial_coords);
 
    // average normal surface
    FTensor::Tensor1<double, 3> t_normal;
    t_normal(i) = t_master_normal(i) - t_slave_normal(i);
    t_normal.normalize();
 
    // get projection operators
    FTensor::Tensor2<double, 3, 3> t_P;
    t_P(i, j) = t_normal(i) * t_normal(j);
    FTensor::Tensor2<double, 3, 3> t_Q;
    t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
 
    constexpr double beta = 0.5;
 
    auto zeta = 1e-8 * face_data_ptr->eleRadius;
    auto alpha1 = ContactOps::alpha_contact_const;
    auto alpha2 = ContactOps::alpha_contact_quadratic;
 
    // this is regularised std::min(d,s)
    auto f_min_gap = [zeta](auto d, auto s) {
      return 0.5 * (d + s - std::sqrt((d - s) * (d - s) + zeta));
    };
    // this derivative of regularised std::min(d,s)
    auto f_diff_min_gap = [zeta](auto d, auto s) {
      return 0.5 * (1 - (d - s) / std::sqrt((d - s) * (d - s) + zeta));
    };
 
    // add penalty on penetration side
    auto f_barrier = [alpha1, alpha2, f_min_gap, f_diff_min_gap](auto g,
                                                                 auto tn) {
      auto d = alpha1 * g;
      auto b1 =
          0.5 * (tn + f_min_gap(d, tn) + f_diff_min_gap(d, tn) * (d - tn));
      auto b2 = alpha2 * f_min_gap(g, 0) * g;
      return b1 - b2;
    };
 
    FTensor::Tensor1<T2, 3> t_gap_vec;
    t_gap_vec(i) = beta * t_spatial_coords(i) +
                   (1 - beta) * t_slave_point_current(i) - t_spatial_coords(i);
    FTensor::Tensor1<typename FTensor::promote<T2, T3>::V, 3> t_traction_vec;
    t_traction_vec(i) =
        -beta * t_master_traction(i) + (beta - 1) * t_slave_traction_current(i);
 
    auto t_gap = t_normal(i) * t_gap_vec(i);
    auto t_tn = t_normal(i) * t_traction_vec(i);
    auto barrier = f_barrier(t_gap, t_tn);
 
    FTensor::Tensor1<typename FTensor::promote<T2, T3>::V, 3> t_traction_bar;
    // add penalty on penetration side
    t_traction_bar(i) = t_normal(i) * barrier;
 
    t_rhs(i) =
 
        t_Q(i, j) * t_master_traction(j) +
 
        t_P(i, j) * (t_master_traction(j) - t_traction_bar(j));
 
    if (debug) {
      auto is_nan_or_inf = [](double value) -> bool {
        return std::isnan(value) || std::isinf(value);
      };
 
      double v = std::complex<double>(t_rhs(i) * t_rhs(i)).real();
      if (is_nan_or_inf(v)) {
        MOFEM_LOG_CHANNEL("SELF");
        MOFEM_LOG("SELF", Sev::error) << "t_rhs " << t_rhs;
        CHK_MOAB_THROW(MOFEM_DATA_INCONSISTENCY, "rhs is nan or inf");
      }
    }
 
  } break;
  case MultiPointRhsType::P:
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "not implemented");
    break;
  }
 
  return t_rhs;
};

multiSlavePoint()

template<typename T1 >

auto EshelbianPlasticity::multiSlavePoint	(	ContactTree::FaceData *	face_data_ptr,
		FTensor::Tensor1< T1, 3 > &	t_spatial_coords
	)

Definition at line 157 of file EshelbianContact.cpp.

  {
 
  return multiPoint(face_data_ptr->unitRay, face_data_ptr->slavePoint,
                    face_data_ptr->slavePointNodes,
                    face_data_ptr->slaveTractionNodes, t_spatial_coords);
};

sort_eigen_vals()

template<int DIM>

auto EshelbianPlasticity::sort_eigen_vals ( FTensor::Tensor1< double, DIM > &  eig, FTensor::Tensor2< double, DIM, DIM > &  eigen_vec  )

inline

Examples

EshelbianOperators.cpp.

Definition at line 159 of file EshelbianOperators.cpp.

                                                                         {
 
  isEq::absMax =
      std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
 
  int i = 0, j = 1, k = 2;
 
  if (is_eq(eig(0), eig(1))) {
    i = 0;
    j = 2;
    k = 1;
  } else if (is_eq(eig(0), eig(2))) {
    i = 0;
    j = 1;
    k = 2;
  } else if (is_eq(eig(1), eig(2))) {
    i = 1;
    j = 0;
    k = 2;
  }
 
  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
      eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
 
      eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
 
      eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
 
  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
 
  {
    FTensor::Index<'i', DIM> i;
    FTensor::Index<'j', DIM> j;
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  }
}

symm_L_tensor()

auto EshelbianPlasticity::symm_L_tensor ( )

inline

Definition at line 52 of file EshelbianAux.hpp.

                            {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'L', size_symm> L;
  FTensor::Dg<double, 3, size_symm> t_L;
  t_L(i, j, L) = 0;
  t_L(0, 0, 0) = 1;
  t_L(1, 1, 3) = 1;
  t_L(2, 2, 5) = 1;
  t_L(1, 0, 1) = 1;
  t_L(2, 0, 2) = 1;
  t_L(2, 1, 4) = 1;
  return t_L;
}

Variable Documentation

size_symm

constexpr static auto EshelbianPlasticity::size_symm = (3 * (3 + 1)) / 2

staticconstexpr

Definition at line 39 of file EshelbianAux.hpp.