EshelbianPlasticity Namespace Reference

Classes
struct	BcDisp
struct	BcRot
struct	CGGUserPolynomialBase
struct	ContactTree
struct	DataAtIntegrationPts
struct	EshelbianCore
struct	FaceRule
struct	FaceUserDataOperatorStabAssembly
struct	FeTractionBc
struct	HMHHencky
struct	HMHPMooneyRivlinWriggersEq63
struct	HMHStVenantKirchhoff
struct	isEq
struct	OpAssembleBasic
struct	OpAssembleFace
struct	OpAssembleVolume
struct	OpCalculateRotationAndSpatialGradient
struct	OpCalculateStrainEnergy
struct	OpConstrainBoundaryHDivLhs_dU
struct	OpConstrainBoundaryHDivRhs
struct	OpConstrainBoundaryL2Lhs_dP
struct	OpConstrainBoundaryL2Lhs_dU
struct	OpConstrainBoundaryL2Rhs
struct	OpDispBc
struct	OpDispBc_dx
struct	OpHMHH
struct	OpJacobian
struct	OpMoveNode
struct	OpPostProcDataStructure
struct	OpRotationBc
struct	OpRotationBc_dx
struct	OpSpatialConsistency_dBubble_domega
struct	OpSpatialConsistency_dP_domega
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	OpSpatialPhysical_du_dx
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	OpTractionBc
struct	OpTreeSearch
struct	PhysicalEquations
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 }
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	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)
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)
auto	diff_deviator (FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
auto	diff_tensor ()
auto	symm_L_tensor ()
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
enum RotSelector	EshelbianCore

Typedef Documentation

BcDispVec

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

Definition at line 303 of file EshelbianPlasticity.hpp.

BcRotVec

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

Definition at line 312 of file EshelbianPlasticity.hpp.

EntData

using EshelbianPlasticity::EntData = typedef EntitiesFieldData::EntData

Definition at line 26 of file EshelbianPlasticity.hpp.

FaceUserDataOperator

using EshelbianPlasticity::FaceUserDataOperator = typedef FaceElementForcesAndSourcesCore::UserDataOperator

Examples

NavierStokesElement.hpp.

Definition at line 29 of file EshelbianPlasticity.hpp.

MatrixPtr

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

Definition at line 23 of file EshelbianPlasticity.hpp.

TractionBcVec

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

Definition at line 323 of file EshelbianPlasticity.hpp.

TractionFreeBc

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

Definition at line 314 of file EshelbianPlasticity.hpp.

UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

Definition at line 27 of file EshelbianPlasticity.hpp.

VectorPtr

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

Definition at line 24 of file EshelbianPlasticity.hpp.

VolUserDataOperator

using EshelbianPlasticity::VolUserDataOperator = typedef VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 28 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

MultiPointRhsType

enum EshelbianPlasticity::MultiPointRhsType

Enumerator
U
P

Definition at line 193 of file EshelbianContact.cpp.

{ U, P };

RotSelector

enum EshelbianPlasticity::RotSelector

Enumerator
SMALL_ROT
MODERATE_ROT
LARGE_ROT

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 20 of file EshelbianPlasticity.hpp.

{ SMALL_ROT, MODERATE_ROT, LARGE_ROT };

StretchSelector

enum EshelbianPlasticity::StretchSelector

Enumerator
LINEAR
LOG

Examples

EshelbianPlasticity.cpp.

Definition at line 21 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 27 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

Examples

EshelbianOperators.cpp.

Definition at line 20 of file EshelbianOperators.cpp.

                                                                    {
  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_tensor()

auto EshelbianPlasticity::diff_tensor ( )

inline

Examples

EshelbianOperators.cpp.

Definition at line 49 of file EshelbianOperators.cpp.

                          {
  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 167 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 119 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 83 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 207 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 36 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 203 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
	)

Definition at line 169 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 143 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
	)

Definition at line 81 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
	)

Definition at line 196 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);
 
    FTensor::Tensor1<double, 3> t_normal;
    t_normal(i) = t_master_normal(i) - t_slave_normal(i);
    t_normal.normalize();
 
    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;
 
    auto f_min_gap = [zeta](auto d, auto s) {
      return 0.5 * (d + s - std::sqrt((d - s) * (d - s) + zeta));
    };
    auto f_diff_min_gap = [zeta](auto d, auto s) {
      return 0.5 * (1 - (d - s) / std::sqrt((d - s) * (d - s) + zeta));
    };
 
    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;
    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 156 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 218 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

Examples

EshelbianOperators.cpp.

Definition at line 60 of file EshelbianOperators.cpp.

                            {
  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

EshelbianCore

enum StretchSelector EshelbianPlasticity::EshelbianCore

Examples

EshelbianPlasticity.cpp.

Definition at line 84 of file EshelbianPlasticity.cpp.

size_symm

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

staticconstexpr

Definition at line 47 of file EshelbianOperators.cpp.