EshelbianPlasticity Namespace Reference

Classes
struct	BcDisp
struct	BcRot
struct	CGGUserPolynomialBase
struct	ContactTree
struct	DataAtIntegrationPts
struct	DynamicRelaxationTimeScale
struct	FaceRule
struct	FaceUserDataOperatorStabAssembly
struct	HMHHencky
struct	HMHNeohookean
struct	HMHPMooneyRivlinWriggersEq63
struct	HMHStVenantKirchhoff
struct	isEq
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	OpGetScale
struct	OpHMHH
struct	OpMoveNode
struct	OpSpatialPhysical
struct	OpSpatialPhysical_du_du
struct	OpTreeSearch
struct	PhysicalEquations
struct	SetIntegrationAtFrontFace
struct	SetIntegrationAtFrontVolume
struct	solve_elastic_setup
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
using	EleOnSide = PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle
using	SideEleOp = EleOnSide::UserDataOperator
typedef std::vector< BcDisp >	BcDispVec
typedef std::vector< BcRot >	BcRotVec
typedef std::vector< Range >	TractionFreeBc
typedef std::vector< TractionBc >	TractionBcVec

Enumerations
enum	SymmetrySelector { SYMMETRIC, NOT_SYMMETRIC, ANIT_SYMMETRIC }
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	is_eq (const double &a, const double &b)
template<int DIM>
auto	get_uniq_nb (double *ptr)
template<typename A , typename B >
auto	sort_eigen_vals (A &eig, B &eigen_vec)
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)

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

Typedef Documentation

BcDispVec

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

Definition at line 424 of file EshelbianPlasticity.hpp.

BcRotVec

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

Definition at line 433 of file EshelbianPlasticity.hpp.

EleOnSide

using EshelbianPlasticity::EleOnSide = typedef PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle

Definition at line 55 of file EshelbianPlasticity.hpp.

EntData

using EshelbianPlasticity::EntData = typedef EntitiesFieldData::EntData

Definition at line 51 of file EshelbianPlasticity.hpp.

FaceUserDataOperator

using EshelbianPlasticity::FaceUserDataOperator = typedef FaceElementForcesAndSourcesCore::UserDataOperator

Examples

NavierStokesElement.hpp.

Definition at line 54 of file EshelbianPlasticity.hpp.

MatrixPtr

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

Definition at line 48 of file EshelbianPlasticity.hpp.

SideEleOp

using EshelbianPlasticity::SideEleOp = typedef EleOnSide::UserDataOperator

Definition at line 56 of file EshelbianPlasticity.hpp.

TractionBcVec

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

Definition at line 444 of file EshelbianPlasticity.hpp.

TractionFreeBc

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

Definition at line 435 of file EshelbianPlasticity.hpp.

UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

Definition at line 52 of file EshelbianPlasticity.hpp.

VectorPtr

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

Definition at line 49 of file EshelbianPlasticity.hpp.

VolUserDataOperator

using EshelbianPlasticity::VolUserDataOperator = typedef VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 53 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

MultiPointRhsType

enum EshelbianPlasticity::MultiPointRhsType

Enumerator
U
P

Definition at line 201 of file EshelbianContact.cpp.

{ U, P };

RotSelector

enum EshelbianPlasticity::RotSelector

Enumerator
SMALL_ROT
MODERATE_ROT
LARGE_ROT
NO_H1_CONFIGURATION

Examples

EshelbianPlasticity.cpp.

Definition at line 45 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 46 of file EshelbianPlasticity.hpp.

{ LINEAR, LOG };

SymmetrySelector

enum EshelbianPlasticity::SymmetrySelector

Enumerator
SYMMETRIC
NOT_SYMMETRIC
ANIT_SYMMETRIC

Examples

EshelbianPlasticity.cpp.

Definition at line 44 of file EshelbianPlasticity.hpp.

{ SYMMETRIC, NOT_SYMMETRIC, ANIT_SYMMETRIC };

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 14 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 69 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 43 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_uniq_nb()

template<int DIM>

auto EshelbianPlasticity::get_uniq_nb ( double *  ptr )

inline

Definition at line 124 of file EshelbianAux.hpp.

                                                        {
  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<double>::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();
  if (EshelbianCore::dynamicRelaxation) {
    ts_time = EshelbianCore::dynamicTime;
    ts_time_step = EshelbianCore::dynamicStep;
  }
 
  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

Definition at line 120 of file EshelbianAux.hpp.

                                                    {
  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 177 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 148 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 86 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 206 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 161 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<typename A , typename B >

auto EshelbianPlasticity::sort_eigen_vals ( A &  eig, B &  eigen_vec  )

inline

Definition at line 135 of file EshelbianAux.hpp.

                                                  {
 
  constexpr int dim = 3;
 
  isEq::absMax =
      std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
  constexpr double eps = std::numeric_limits<double>::epsilon();
  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
 
  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(dim, i);
    FTENSOR_INDEX(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 54 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 41 of file EshelbianAux.hpp.