v0.14.0
Classes | Typedefs | Enumerations | Functions | Variables
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< BcDispBcDispVec
 
typedef std::vector< BcRotBcRotVec
 
typedef std::vector< RangeTractionFreeBc
 
typedef std::vector< TractionBcTractionBcVec
 

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

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

Definition at line 323 of file EshelbianPlasticity.hpp.

◆ 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

Definition at line 28 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

◆ MultiPointRhsType

Enumerator

Definition at line 193 of file EshelbianContact.cpp.

193 { U, P };

◆ RotSelector

Enumerator
SMALL_ROT 
MODERATE_ROT 
LARGE_ROT 
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 20 of file EshelbianPlasticity.hpp.

◆ StretchSelector

Enumerator
LINEAR 
LOG 
Examples
EshelbianPlasticity.cpp.

Definition at line 21 of file EshelbianPlasticity.hpp.

21 { 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
ppolynomial order
Nshape functions
diffNdirevatives of shape functions
l2_basebase functions for l2 space
diff_l2_basedirevatives base functions for l2 space
phireturned base functions
gdimnumber of integration points
Returns
MoFEMErrorCode
Examples
EshelbianPlasticity.cpp.

Definition at line 20 of file CGGTonsorialBubbleBase.cpp.

23  {
25 
26  // FIXME: In this implementation symmetry of mix derivatives is not exploited
27 
28  if (p < 1)
30 
36 
37  Tensor1<double, 3> t_diff_n[4];
38  {
39  Tensor1<PackPtr<const double *, 3>, 3> t_diff_n_tmp(&diffN[0], &diffN[1],
40  &diffN[2]);
41  for (int ii = 0; ii != 4; ++ii) {
42  t_diff_n[ii](i) = t_diff_n_tmp(i);
43  ++t_diff_n_tmp;
44  }
45  }
46 
47  Tensor1<double, 3> t_diff_ksi[3];
48  for (int ii = 0; ii != 3; ++ii)
49  t_diff_ksi[ii](i) = t_diff_n[ii + 1](i) - t_diff_n[0](i);
50 
51  int lp = p >= 2 ? p - 2 + 1 : 0;
52  VectorDouble l[3] = {VectorDouble(lp + 1), VectorDouble(lp + 1),
53  VectorDouble(lp + 1)};
54  MatrixDouble diff_l[3] = {MatrixDouble(3, lp + 1), MatrixDouble(3, lp + 1),
55  MatrixDouble(3, lp + 1)};
56  MatrixDouble diff2_l[3] = {MatrixDouble(9, lp + 1), MatrixDouble(9, lp + 1),
57  MatrixDouble(9, lp + 1)};
58 
59  for (int ii = 0; ii != 3; ++ii)
60  diff2_l[ii].clear();
61 
62  for (int gg = 0; gg != gdim; ++gg) {
63 
64  const int node_shift = gg * 4;
65 
66  for (int ii = 0; ii != 3; ++ii) {
67 
68  auto &t_diff_ksi_ii = t_diff_ksi[ii];
69  auto &l_ii = l[ii];
70  auto &diff_l_ii = diff_l[ii];
71  auto &diff2_l_ii = diff2_l[ii];
72 
73  double ksi_ii = N[node_shift + ii + 1] - N[node_shift + 0];
74 
75  CHKERR Legendre_polynomials(lp, ksi_ii, &t_diff_ksi_ii(0),
76  &*l_ii.data().begin(),
77  &*diff_l_ii.data().begin(), 3);
78 
79  for (int l = 1; l < lp; ++l) {
80  const double a = ((2 * (double)l + 1) / ((double)l + 1));
81  const double b = ((double)l / ((double)l + 1));
82  for (int d0 = 0; d0 != 3; ++d0)
83  for (int d1 = 0; d1 != 3; ++d1) {
84  const int r = 3 * d0 + d1;
85  diff2_l_ii(r, l + 1) = a * (t_diff_ksi_ii(d0) * diff_l_ii(d1, l) +
86  t_diff_ksi_ii(d1) * diff_l_ii(d0, l) +
87  ksi_ii * diff2_l_ii(r, l)) -
88  b * diff2_l_ii(r, l - 1);
89  }
90  }
91  }
92 
93  const double n[] = {N[node_shift + 0], N[node_shift + 1], N[node_shift + 2],
94  N[node_shift + 3]};
95 
97  Tensor3<double, 3, 3, 3> t_bk_diff;
98  const int tab[4][4] = {
99  {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
100  t_bk(i, j) = 0;
101  t_bk_diff(i, j, k) = 0;
102  for (int ii = 0; ii != 3; ++ii) {
103  const int i0 = tab[ii][0];
104  const int i1 = tab[ii][1];
105  const int i2 = tab[ii][2];
106  const int i3 = tab[ii][3];
107  auto &t_diff_n_i0 = t_diff_n[i0];
108  auto &t_diff_n_i1 = t_diff_n[i1];
109  auto &t_diff_n_i2 = t_diff_n[i2];
110  auto &t_diff_n_i3 = t_diff_n[i3];
112  t_k(i, j) = t_diff_n_i3(i) * t_diff_n_i3(j);
113  const double b = n[i0] * n[i1] * n[i2];
114  t_bk(i, j) += b * t_k(i, j);
115  Tensor1<double, 3> t_diff_b;
116  t_diff_b(i) = t_diff_n_i0(i) * n[i1] * n[i2] +
117  t_diff_n_i1(i) * n[i0] * n[i2] +
118  t_diff_n_i2(i) * n[i0] * n[i1];
119  t_bk_diff(i, j, k) += t_k(i, j) * t_diff_b(k);
120  }
121 
122  int zz = 0;
123  for (int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
124 
125  for (int ii = 0; ii <= o; ++ii)
126  for (int jj = 0; (ii + jj) <= o; ++jj) {
127 
128  const int kk = o - ii - jj;
129 
130  auto get_diff_l = [&](const int y, const int i) {
131  return Tensor1<double, 3>(diff_l[y](0, i), diff_l[y](1, i),
132  diff_l[y](2, i));
133  };
134  auto get_diff2_l = [&](const int y, const int i) {
135  return Tensor2<double, 3, 3>(
136  diff2_l[y](0, i), diff2_l[y](1, i), diff2_l[y](2, i),
137  diff2_l[y](3, i), diff2_l[y](4, i), diff2_l[y](5, i),
138  diff2_l[y](6, i), diff2_l[y](7, i), diff2_l[y](8, i));
139  };
140 
141  auto l_i = l[0][ii];
142  auto t_diff_i = get_diff_l(0, ii);
143  auto t_diff2_i = get_diff2_l(0, ii);
144  auto l_j = l[1][jj];
145  auto t_diff_j = get_diff_l(1, jj);
146  auto t_diff2_j = get_diff2_l(1, jj);
147  auto l_k = l[2][kk];
148  auto t_diff_k = get_diff_l(2, kk);
149  auto t_diff2_k = get_diff2_l(2, kk);
150 
151  Tensor1<double, 3> t_diff_l2;
152  t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
153  t_diff_k(i) * l_i * l_j;
154  Tensor2<double, 3, 3> t_diff2_l2;
155  t_diff2_l2(i, j) =
156  t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
157  t_diff_i(i) * l_j * t_diff_k(j) +
158 
159  t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
160  t_diff_j(i) * l_i * t_diff_k(j) +
161 
162  t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
163  t_diff_k(i) * l_i * t_diff_j(j);
164 
165  for (int dd = 0; dd != 3; ++dd) {
166 
167  Tensor2<double, 3, 3> t_axial_diff;
168  t_axial_diff(i, j) = 0;
169  for (int mm = 0; mm != 3; ++mm)
170  t_axial_diff(dd, mm) = t_diff_l2(mm);
171 
172  Tensor3<double, 3, 3, 3> t_A_diff;
173  t_A_diff(i, j, k) = levi_civita(i, j, m) * t_axial_diff(m, k);
174  Tensor2<double, 3, 3> t_curl_A;
175  t_curl_A(i, j) = levi_civita(j, m, f) * t_A_diff(i, f, m);
176  Tensor3<double, 3, 3, 3> t_curl_A_bK_diff;
177  t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
178 
179  Tensor3<double, 3, 3, 3> t_axial_diff2;
180  t_axial_diff2(i, j, k) = 0;
181  for (int mm = 0; mm != 3; ++mm)
182  for (int nn = 0; nn != 3; ++nn)
183  t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
184  Tensor4<double, 3, 3, 3, 3> t_A_diff2;
185  t_A_diff2(i, j, k, f) =
186  levi_civita(i, j, m) * t_axial_diff2(m, k, f);
187  Tensor3<double, 3, 3, 3> t_curl_A_diff2;
188  t_curl_A_diff2(i, j, k) =
189  levi_civita(j, m, f) * t_A_diff2(i, f, m, k);
190  Tensor3<double, 3, 3, 3> t_curl_A_diff2_bK;
191  t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
192 
193  t_phi(i, j) = levi_civita(j, m, f) * (t_curl_A_bK_diff(i, f, m) +
194  t_curl_A_diff2_bK(i, f, m));
195 
196  ++t_phi;
197  ++zz;
198  }
199  }
200  }
201  if (zz != NBVOLUMETET_CCG_BUBBLE(p))
202  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
203  "Wrong number of base functions %d != %d", zz,
205  }
206 
208 }

◆ checkSdf()

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

Definition at line 27 of file EshelbianContact.cpp.

27  {
28  for (auto &m_sdf : sdf_map_range) {
29  if (m_sdf.second.find(fe_ent) != m_sdf.second.end())
30  return m_sdf.first;
31  }
32  return -1;
33 }

◆ 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.

20  {
25 
26  FTensor::Ddg<double, 3, 3> t_diff_deviator;
27  t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
28 
29  constexpr double third = boost::math::constants::third<double>();
30 
31  t_diff_deviator(0, 0, 0, 0) -= third;
32  t_diff_deviator(0, 0, 1, 1) -= third;
33 
34  t_diff_deviator(1, 1, 0, 0) -= third;
35  t_diff_deviator(1, 1, 1, 1) -= third;
36 
37  t_diff_deviator(2, 2, 0, 0) -= third;
38  t_diff_deviator(2, 2, 1, 1) -= third;
39 
40  t_diff_deviator(0, 0, 2, 2) -= third;
41  t_diff_deviator(1, 1, 2, 2) -= third;
42  t_diff_deviator(2, 2, 2, 2) -= third;
43 
44  return t_diff_deviator;
45 }

◆ diff_tensor()

auto EshelbianPlasticity::diff_tensor ( )
inline
Examples
EshelbianOperators.cpp.

Definition at line 49 of file EshelbianOperators.cpp.

49  {
56  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
57  return t_diff;
58 };

◆ 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.

168  {
169 
170  using D = typename TensorTypeExtractor<T>::Type;
171 
174 
175  constexpr double eps = 1e-10;
176  for (int l = 0; l != 3; ++l) {
178  t_omega_c(i) = t_omega(i);
179  t_omega_c(l) += std::complex<double>(0, eps);
180  auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
181  for (int i = 0; i != 3; ++i) {
182  for (int j = 0; j != 3; ++j) {
183  for (int k = 0; k != 3; ++k) {
184  t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
185  }
186  }
187  }
188  }
189 
190  return t_diff2_R;
191 }

◆ 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.

120  {
121 
122  using D = typename TensorTypeExtractor<T>::Type;
123 
128 
130 
131  const auto angle =
132  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
133  if (rotSelector == SMALL_ROT) {
134  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
135  return t_diff_R;
136  }
137 
138  const auto ss = sin(angle);
139  const auto a = ss / angle;
140  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
141 
143  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
144 
145  const auto angle2 = angle * angle;
146  const auto cc = cos(angle);
147  const auto diff_a = (angle * cc - ss) / angle2;
148  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
149  if (rotSelector == MODERATE_ROT)
150  return t_diff_R;
151 
152  const auto ss_2 = sin(angle / 2.);
153  const auto cc_2 = cos(angle / 2.);
154  const auto b = 2. * ss_2 * ss_2 / angle2;
155  t_diff_R(i, j, k) +=
156  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
157  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
158  const auto diff_b =
159  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
160  t_diff_R(i, j, k) +=
161  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
162 
163  return t_diff_R;
164 }

◆ 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.

84  {
85 
86  using D = typename TensorTypeExtractor<T>::Type;
87 
92 
93  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
94  t_R(i, j) = t_kd(i, j);
95 
97  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
98 
99  if (rotSelector == SMALL_ROT) {
100  t_R(i, j) += t_Omega(i, j);
101  return t_R;
102  }
103 
104  const auto angle =
105  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
106  const auto a = sin(angle) / angle;
107  t_R(i, j) += a * t_Omega(i, j);
108  if (rotSelector == MODERATE_ROT)
109  return t_R;
110 
111  const auto ss_2 = sin(angle / 2.) / angle;
112  const auto b = 2. * ss_2 * ss_2;
113  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
114 
115  return t_R;
116 }

◆ 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.

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

◆ 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.

37  {
38 
39  auto nb_gauss_pts = op_ptr->getGaussPts().size2();
40 
41  auto ts_time = op_ptr->getTStime();
42  auto ts_time_step = op_ptr->getTStimeStep();
43 
44  auto m_spatial_coords = ContactOps::get_spatial_coords(
45  op_ptr->getFTensor1CoordsAtGaussPts(),
46  getFTensor1FromMat<3>(contact_disp), nb_gauss_pts);
47  auto m_normals_at_pts = ContactOps::get_normalize_normals(
48  op_ptr->getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
50  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
51  block_id);
53  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
54  block_id);
55 
56 #ifndef NDEBUG
57  if (v_sdf.size() != nb_gauss_pts)
59  "Wrong number of integration pts");
60  if (m_grad_sdf.size2() != nb_gauss_pts)
62  "Wrong number of integration pts");
63  if (m_grad_sdf.size1() != 3)
64  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Should be size of 3");
65 #endif // NDEBUG
66 
67  if (eval_hessian) {
69  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
70  block_id);
71  return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
72  m_hess_sdf);
73  } else {
74 
75  return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
76  MatrixDouble(6, nb_gauss_pts, 0.));
77  }
78 };

◆ is_eq()

auto EshelbianPlasticity::is_eq ( const double a,
const double b 
)
inline
Examples
EshelbianOperators.cpp.

Definition at line 203 of file EshelbianOperators.cpp.

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

◆ 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.

170  {
171 
172  auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
173  multiSlavePoint(face_data_ptr, t_spatial_coords);
174  auto [t_master_point_current, t_master_traction_current, t_master_normal] =
175  multiMasterPoint(face_data_ptr, t_spatial_coords);
176 
178 
180  t_normal(i) = t_master_normal(i) - t_slave_normal(i);
181  t_normal.normalize();
182 
183  auto gap = t_normal(i) * (t_slave_point_current(i) - t_spatial_coords(i));
184  auto tn_master = t_master_traction_current(i) * t_normal(i);
185  auto tn_slave = t_slave_traction_current(i) * t_normal(i);
186  auto tn = std::max(-tn_master, tn_slave);
187 
188  return std::make_tuple(gap, tn_master, tn_slave,
189  ContactOps::constrain(gap, tn),
190  t_master_traction_current, t_slave_traction_current);
191 };

◆ 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.

148  {
149 
150  return multiPoint(face_data_ptr->unitRay, face_data_ptr->masterPoint,
151  face_data_ptr->masterPointNodes,
152  face_data_ptr->masterTractionNodes, t_spatial_coords);
153 };

◆ 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.

88  {
90 
91  auto t_unit_ray = getFTensor1FromPtr<3>(unit_ray.data());
92  auto t_point = getFTensor1FromPtr<3>(point.data());
93 
94  auto get_normal = [](auto &ele_coords) {
96  Tools::getTriNormal(ele_coords.data(), &t_normal(0));
97  return t_normal;
98  };
99 
100  auto t_normal = get_normal(elem_point_nodes);
101  t_normal(i) /= t_normal.l2();
102 
103  auto sn = t_normal(i) * t_point(i);
104  auto nm = t_normal(i) * t_spatial_coords(i);
105  auto nr = t_normal(i) * t_unit_ray(i);
106 
107  auto gamma = (sn - nm) / nr;
108 
109  FTensor::Tensor1<T1, 3> t_point_current;
110  t_point_current(i) = t_spatial_coords(i) + gamma * t_unit_ray(i);
111 
112  auto get_local_point_shape_functions = [&](auto &&t_elem_coords,
113  auto &t_point) {
114  std::array<T1, 2> loc_coords;
116  Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
117  &t_elem_coords(0, 0), &t_point(0), 1, loc_coords.data()),
118  "get local coords");
119  return FTensor::Tensor1<T1, 3>{N_MBTRI0(loc_coords[0], loc_coords[1]),
120  N_MBTRI1(loc_coords[0], loc_coords[1]),
121  N_MBTRI2(loc_coords[0], loc_coords[1])};
122  };
123 
124  auto eval_position = [&](auto &&t_field, auto &t_shape_fun) {
127  FTensor::Tensor1<T1, 3> t_point_out;
128  t_point_out(i) = t_shape_fun(j) * t_field(j, i);
129  return t_point_out;
130  };
131 
132  auto t_shape_fun = get_local_point_shape_functions(
133  getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_point_current);
134  auto t_slave_point_updated = eval_position(
135  getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_shape_fun);
136  auto t_traction_updated = eval_position(
137  getFTensor2FromPtr<3, 3>(elem_traction_nodes.data()), t_shape_fun);
138 
139  return std::make_tuple(t_slave_point_updated, t_traction_updated, t_normal);
140 };

◆ 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.

203  {
206 
208  t_u(i) = t_spatial_coords(i) - t_coords(i);
209 
211 
212  switch (type) {
213  case MultiPointRhsType::U: {
214 
215  auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
216  multiSlavePoint(face_data_ptr, t_spatial_coords);
217  auto [t_master_point_current, t_master_traction_current, t_master_normal] =
218  multiMasterPoint(face_data_ptr, t_spatial_coords);
219 
221  t_normal(i) = t_master_normal(i) - t_slave_normal(i);
222  t_normal.normalize();
223 
225  t_P(i, j) = t_normal(i) * t_normal(j);
227  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
228 
229  constexpr double beta = 0.5;
230 
231  auto zeta = 1e-8 * face_data_ptr->eleRadius;
232  auto alpha1 = ContactOps::alpha_contact_const;
234 
235  auto f_min_gap = [zeta](auto d, auto s) {
236  return 0.5 * (d + s - std::sqrt((d - s) * (d - s) + zeta));
237  };
238  auto f_diff_min_gap = [zeta](auto d, auto s) {
239  return 0.5 * (1 - (d - s) / std::sqrt((d - s) * (d - s) + zeta));
240  };
241 
242  auto f_barrier = [alpha1, alpha2, f_min_gap, f_diff_min_gap](auto g,
243  auto tn) {
244  auto d = alpha1 * g;
245  auto b1 =
246  0.5 * (tn + f_min_gap(d, tn) + f_diff_min_gap(d, tn) * (d - tn));
247  auto b2 = alpha2 * f_min_gap(g, 0) * g;
248  return b1 - b2;
249  };
250 
251  FTensor::Tensor1<T2, 3> t_gap_vec;
252  t_gap_vec(i) = beta * t_spatial_coords(i) +
253  (1 - beta) * t_slave_point_current(i) - t_spatial_coords(i);
255  t_traction_vec(i) =
256  -beta * t_master_traction(i) + (beta - 1) * t_slave_traction_current(i);
257 
258  auto t_gap = t_normal(i) * t_gap_vec(i);
259  auto t_tn = t_normal(i) * t_traction_vec(i);
260  auto barrier = f_barrier(t_gap, t_tn);
261 
263  t_traction_bar(i) = t_normal(i) * barrier;
264 
265  t_rhs(i) =
266 
267  t_Q(i, j) * t_master_traction(j) +
268 
269  t_P(i, j) * (t_master_traction(j) - t_traction_bar(j));
270 
271  if (debug) {
272  auto is_nan_or_inf = [](double value) -> bool {
273  return std::isnan(value) || std::isinf(value);
274  };
275 
276  double v = std::complex<double>(t_rhs(i) * t_rhs(i)).real();
277  if (is_nan_or_inf(v)) {
278  MOFEM_LOG_CHANNEL("SELF");
279  MOFEM_LOG("SELF", Sev::error) << "t_rhs " << t_rhs;
280  CHK_MOAB_THROW(MOFEM_DATA_INCONSISTENCY, "rhs is nan or inf");
281  }
282  }
283 
284  } break;
286  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "not implemented");
287  break;
288  }
289 
290  return t_rhs;
291 };

◆ 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.

161  {
162 
163  return multiPoint(face_data_ptr->unitRay, face_data_ptr->slavePoint,
164  face_data_ptr->slavePointNodes,
165  face_data_ptr->slaveTractionNodes, t_spatial_coords);
166 };

◆ 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.

219  {
220 
221  isEq::absMax =
222  std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
223 
224  int i = 0, j = 1, k = 2;
225 
226  if (is_eq(eig(0), eig(1))) {
227  i = 0;
228  j = 2;
229  k = 1;
230  } else if (is_eq(eig(0), eig(2))) {
231  i = 0;
232  j = 1;
233  k = 2;
234  } else if (is_eq(eig(1), eig(2))) {
235  i = 1;
236  j = 0;
237  k = 2;
238  }
239 
240  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
241  eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
242 
243  eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
244 
245  eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
246 
247  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
248 
249  {
252  eig(i) = eig_c(i);
253  eigen_vec(i, j) = eigen_vec_c(i, j);
254  }
255 }

◆ symm_L_tensor()

auto EshelbianPlasticity::symm_L_tensor ( )
inline
Examples
EshelbianOperators.cpp.

Definition at line 60 of file EshelbianOperators.cpp.

60  {
65  t_L(i, j, L) = 0;
66  t_L(0, 0, 0) = 1;
67  t_L(1, 1, 3) = 1;
68  t_L(2, 2, 5) = 1;
69  t_L(1, 0, 1) = 1;
70  t_L(2, 0, 2) = 1;
71  t_L(2, 1, 4) = 1;
72  return t_L;
73 }

Variable Documentation

◆ 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
ContactOps::get_spatial_coords
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
Definition: ContactOps.hpp:407
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
EshelbianPlasticity::is_eq
auto is_eq(const double &a, const double &b)
Definition: EshelbianOperators.cpp:203
g
constexpr double g
Definition: shallow_wave.cpp:63
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:21
FTensor::Tensor1< double, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
ContactOps::hess_surface_distance_function
MatrixDouble hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:332
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:574
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
a
constexpr double a
Definition: approx_sphere.cpp:30
ContactOps::get_normalize_normals
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:424
double
convert.type
type
Definition: convert.py:64
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
EshelbianPlasticity::multiPoint
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)
Definition: EshelbianContact.cpp:81
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
EshelbianPlasticity::multiMasterPoint
auto multiMasterPoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Definition: EshelbianContact.cpp:143
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
ContactOps::alpha_contact_const
double alpha_contact_const
Definition: EshelbianContact.hpp:20
FTensor::Index< 'i', 3 >
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
convert.n
n
Definition: convert.py:82
N_MBTRI0
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:46
N
const int N
Definition: speed_test.cpp:3
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Dg
Definition: Dg_value.hpp:9
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:21
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
N_MBTRI1
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:47
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
FTensor::Ddg< double, 3, 3 >
EshelbianPlasticity::get_diff_rotation_form_vector
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:119
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
ContactOps::surface_distance_function
VectorDouble surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:216
EshelbianPlasticity::multiSlavePoint
auto multiSlavePoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Definition: EshelbianContact.cpp:156
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
ContactOps::grad_surface_distance_function
MatrixDouble grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:275
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ContactOps::alpha_contact_quadratic
double alpha_contact_quadratic
Definition: EshelbianContact.hpp:21
EshelbianPlasticity::U
@ U
Definition: EshelbianContact.cpp:193
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
N_MBTRI2
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:48