v0.10.0
Classes | Typedefs | Enumerations | Functions | Variables
EshelbianPlasticity Namespace Reference

Classes

struct  BcDisp
 
struct  BcRot
 
struct  CGGUserPolynomialBase
 
struct  DataAtIntegrationPts
 
struct  EpElement
 
struct  EpElement< BasicMethod >
 
struct  EpElement< FEMethod >
 
struct  EpElement< FeTractionBc >
 
struct  EpElementBase
 
struct  EpFEMethod
 
struct  EshelbianCore
 
struct  FaceRule
 
struct  FeTractionBc
 
struct  HMHPMooneyRivlinWriggersEq63
 
struct  HMHStVenantKirchhoff
 
struct  OpAssembleBasic
 
struct  OpAssembleFace
 
struct  OpAssembleVolume
 
struct  OpCalculateRotationAndSpatialGradient
 
struct  OpCalculateStrainEnergy
 
struct  OpDispBc
 
struct  OpDispBc_dx
 
struct  OpHMHH
 
struct  OpJacobian
 
struct  OpL2Transform
 
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  OpSpatialPreconditionMass
 
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  OpSpatialSchurBegin
 
struct  OpSpatialSchurEnd
 
struct  OpTractionBc
 
struct  PhysicalEquations
 
struct  TensorTypeExtractor
 
struct  TensorTypeExtractor< FTensor::PackPtr< T, I > >
 
struct  TractionBc
 
struct  VolRule
 Set integration rule on element. More...
 

Typedefs

typedef boost::shared_ptr< MatrixDouble > MatrixPtr
 
typedef boost::shared_ptr< VectorDouble > VectorPtr
 
using EntData = DataForcesAndSourcesCore::EntData
 
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator
 
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
 
using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator
 
typedef std::vector< BcDispBcDispVec
 
typedef std::vector< BcRotBcRotVec
 
typedef std::vector< Range > TractionFreeBc
 
typedef std::vector< TractionBcTractionBcVec
 

Enumerations

enum  EshelbianInterfaces { CORE_ESHELBIAN_INTERFACE = 1 << 0 }
 

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...
 
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat (MatrixDouble &m)
 
template<typename T >
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector (FTensor::Tensor1< T, 3 > &t_omega)
 
template<typename T >
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector (FTensor::Tensor1< T, 3 > &t_omega)
 
template<typename T >
static MoFEMErrorCode tensor_exponent (FTensor::Tensor2_symmetric< T, 3 > &t_log_s_u, FTensor::Tensor2_symmetric< T, 3 > &t_s_u, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff_s_u)
 
static MoFEMErrorCode matrix_exponent_derivative (MatrixDouble &log_u, MatrixDouble &diff_u, std::array< MatrixDouble, 6 > &diff2_u)
 
static MoFEMErrorCode matrix_rotation_derivative (MatrixDouble &omega, std::array< MatrixDouble, 3 > &diff2_omega)
 
template<typename T >
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > calculate_viscous_stress (FTensor::Tensor2< T, 3, 3 > &t)
 
template<typename T >
static FTensor::Tensor4< typename TensorTypeExtractor< T >::Type, 3, 3, 3, 3 > calculate_diff_viscous_stress (FTensor::Tensor2< T, 3, 3 > &t, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff)
 

Variables

static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
 

Typedef Documentation

◆ BcDispVec

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

Definition at line 589 of file EshelbianPlasticity.hpp.

◆ BcRotVec

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

Definition at line 598 of file EshelbianPlasticity.hpp.

◆ EntData

using EshelbianPlasticity::EntData = typedef DataForcesAndSourcesCore::EntData

Definition at line 34 of file EshelbianPlasticity.hpp.

◆ FaceUserDataOperator

Definition at line 37 of file EshelbianPlasticity.hpp.

◆ MatrixPtr

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

Definition at line 20 of file EshelbianPlasticity.hpp.

◆ TractionBcVec

Definition at line 609 of file EshelbianPlasticity.hpp.

◆ TractionFreeBc

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

Definition at line 600 of file EshelbianPlasticity.hpp.

◆ UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

Definition at line 35 of file EshelbianPlasticity.hpp.

◆ VectorPtr

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

Definition at line 21 of file EshelbianPlasticity.hpp.

◆ VolUserDataOperator

Definition at line 36 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

◆ EshelbianInterfaces

Enumerator
CORE_ESHELBIAN_INTERFACE 

Definition at line 15 of file EshelbianPlasticity.hpp.

Function Documentation

◆ calculate_diff_viscous_stress()

template<typename T >
static FTensor::Tensor4<typename TensorTypeExtractor<T>::Type, 3, 3, 3, 3> EshelbianPlasticity::calculate_diff_viscous_stress ( FTensor::Tensor2< T, 3, 3 > &  t,
FTensor::Tensor4< T, 3, 3, 3, 3 > &  t_diff 
)
static
Examples
EshelbianOperators.cpp.

Definition at line 379 of file EshelbianOperators.cpp.

380  {
381 
387 
388  auto get_dot_fun = [i, j, k, m, n](auto &t, auto &t_diff) {
390  t_diff_c;
391  t_diff_c(i, j, m, n) = t_diff(i, j, m, n);
392  return t_diff_c;
393  };
394 
395  auto get_trace = [i, j](auto &t_diff) {
397  t_diff_tr(i, j) = 0;
398  for (auto ii = 0; ii != 3; ++ii)
399  for (auto mm = 0; mm != 3; ++mm)
400  for (auto nn = 0; nn != 3; ++nn)
401  t_diff_tr(mm, nn) += t_diff(ii, ii, mm, nn) / 3.;
402  return t_diff_tr;
403  };
404 
405  auto get_deviator = [i, j, m, n, get_trace](auto &&t_diff) {
406  auto t_diff_tr = get_trace(t_diff);
408  t_diff_dev;
409  t_diff_dev(i, j, m, n) = t_diff(i, j, m, n);
410  for (auto mm = 0; mm != 3; ++mm)
411  for (auto nn = 0; nn != 3; ++nn)
412  for (auto ii = 0; ii != 3; ++ii)
413  t_diff_dev(ii, ii, mm, nn) -= t_diff_tr(mm, nn);
414  return t_diff_dev;
415  };
416 
417  return get_deviator(get_dot_fun(t, t_diff));
418 }
static Index< 'n', 3 > n
static Index< 'm', 3 > m
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ calculate_viscous_stress()

template<typename T >
static FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3, 3> EshelbianPlasticity::calculate_viscous_stress ( FTensor::Tensor2< T, 3, 3 > &  t)
static
Examples
EshelbianOperators.cpp.

Definition at line 346 of file EshelbianOperators.cpp.

346  {
347 
351 
352  auto get_dot_fun = [i, j, k](auto &t) {
354  t_c(i, j) = t(i, j);
355  return t_c;
356  };
357 
358  auto get_trace = [](auto &t) {
359  typename TensorTypeExtractor<T>::Type tr = 0;
360  for (auto ii = 0; ii != 3; ++ii)
361  tr += t(ii, ii) / 3.;
362  return tr;
363  };
364 
365  auto get_deviator = [i, j, get_trace](auto &&t) {
366  const auto tr = get_trace(t);
368  t_dev(i, j) = t(i, j);
369  for (auto ii = 0; ii != 3; ++ii)
370  t_dev(ii, ii) -= tr;
371  return t_dev;
372  };
373 
374  return get_deviator(get_dot_fun(t));
375 }
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ 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 [14]

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 
35  Index<'f', 3> f;
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 
94  const double n[] = {N[node_shift + 0], N[node_shift + 1], N[node_shift + 2],
95  N[node_shift + 3]};
96 
98  Tensor3<double, 3, 3, 3> t_bk_diff;
99  const int tab[4][4] = {
100  {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
101  t_bk(i, j) = 0;
102  t_bk_diff(i, j, k) = 0;
103  for (int ii = 0; ii != 3; ++ii) {
104  const int i0 = tab[ii][0];
105  const int i1 = tab[ii][1];
106  const int i2 = tab[ii][2];
107  const int i3 = tab[ii][3];
108  auto &t_diff_n_i0 = t_diff_n[i0];
109  auto &t_diff_n_i1 = t_diff_n[i1];
110  auto &t_diff_n_i2 = t_diff_n[i2];
111  auto &t_diff_n_i3 = t_diff_n[i3];
113  t_k(i, j) = t_diff_n_i3(i) * t_diff_n_i3(j);
114  const double b = n[i0] * n[i1] * n[i2];
115  t_bk(i, j) += b * t_k(i, j);
116  Tensor1<double, 3> t_diff_b;
117  t_diff_b(i) = t_diff_n_i0(i) * n[i1] * n[i2] +
118  t_diff_n_i1(i) * n[i0] * n[i2] +
119  t_diff_n_i2(i) * n[i0] * n[i1];
120  t_bk_diff(i, j, k) += t_k(i, j) * t_diff_b(k);
121  }
122 
123  int zz = 0;
124  for (int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
125 
126  for (int ii = 0; ii <= o; ++ii)
127  for (int jj = 0; (ii + jj) <= o; ++jj) {
128 
129  const int kk = o - ii - jj;
130 
131  auto get_diff_l = [&](const int y, const int i) {
132  return Tensor1<double, 3>(diff_l[y](0, i), diff_l[y](1, i),
133  diff_l[y](2, i));
134  };
135  auto get_diff2_l = [&](const int y, const int i) {
136  return Tensor2<double, 3, 3>(
137  diff2_l[y](0, i), diff2_l[y](1, i), diff2_l[y](2, i),
138  diff2_l[y](3, i), diff2_l[y](4, i), diff2_l[y](5, i),
139  diff2_l[y](6, i), diff2_l[y](7, i), diff2_l[y](8, i));
140  };
141 
142  auto l_i = l[0][ii];
143  auto t_diff_i = get_diff_l(0, ii);
144  auto t_diff2_i = get_diff2_l(0, ii);
145  auto l_j = l[1][jj];
146  auto t_diff_j = get_diff_l(1, jj);
147  auto t_diff2_j = get_diff2_l(1, jj);
148  auto l_k = l[2][kk];
149  auto t_diff_k = get_diff_l(2, kk);
150  auto t_diff2_k = get_diff2_l(2, kk);
151 
152  Tensor1<double, 3> t_diff_l2;
153  t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
154  t_diff_k(i) * l_i * l_j;
155  Tensor2<double, 3, 3> t_diff2_l2;
156  t_diff2_l2(i, j) =
157  t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
158  t_diff_i(i) * l_j * t_diff_k(j) +
159 
160  t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
161  t_diff_j(i) * l_i * t_diff_k(j) +
162 
163  t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
164  t_diff_k(i) * l_i * t_diff_j(j);
165 
166  for (int dd = 0; dd != 3; ++dd) {
167 
168  Tensor2<double, 3, 3> t_axial_diff;
169  t_axial_diff(i, j) = 0;
170  for (int mm = 0; mm != 3; ++mm)
171  t_axial_diff(dd, mm) = t_diff_l2(mm);
172 
173  Tensor3<double, 3, 3, 3> t_A_diff;
174  t_A_diff(i, j, k) = levi_civita(i, j, m) * t_axial_diff(m, k);
175  Tensor2<double, 3, 3> t_curl_A;
176  t_curl_A(i, j) = levi_civita(j, m, f) * t_A_diff(i, f, m);
177  Tensor3<double, 3, 3, 3> t_curl_A_bK_diff;
178  t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
179 
180  Tensor3<double, 3, 3, 3> t_axial_diff2;
181  t_axial_diff2(i, j, k) = 0;
182  for (int mm = 0; mm != 3; ++mm)
183  for (int nn = 0; nn != 3; ++nn)
184  t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
185  Tensor4<double, 3, 3, 3, 3> t_A_diff2;
186  t_A_diff2(i, j, k, f) =
187  levi_civita(i, j, m) * t_axial_diff2(m, k, f);
188  Tensor3<double, 3, 3, 3> t_curl_A_diff2;
189  t_curl_A_diff2(i, j, k) =
190  levi_civita(j, m, f) * t_A_diff2(i, f, m, k);
191  Tensor3<double, 3, 3, 3> t_curl_A_diff2_bK;
192  t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
193 
194  t_phi(i, j) = levi_civita(j, m, f) * (t_curl_A_bK_diff(i, f, m) +
195  t_curl_A_diff2_bK(i, f, m));
196 
197  ++t_phi;
198  ++zz;
199  }
200  }
201  }
202  if(zz != NBVOLUMETET_CCG_BUBBLE(p))
203  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
204  "Wrong number of base functions %d != %d", zz,
206 
207  }
208 
210 }
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
static Index< 'p', 3 > p
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
static Index< 'n', 3 > n
static Index< 'm', 3 > m
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
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'o', 3 > o
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
static Index< 'k', 3 > k
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
const int N
Definition: speed_test.cpp:3
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

◆ get_diff_rotation_form_vector()

template<typename T >
static FTensor::Tensor3<typename TensorTypeExtractor<T>::Type, 3, 3, 3> EshelbianPlasticity::get_diff_rotation_form_vector ( FTensor::Tensor1< T, 3 > &  t_omega)
static
Examples
EshelbianOperators.cpp.

Definition at line 60 of file EshelbianOperators.cpp.

61  {
62 
67 
68  const typename TensorTypeExtractor<T>::Type angle =
69  sqrt(t_omega(i) * t_omega(i));
70 
71  constexpr double tol = 1e-18;
72  if (std::abs(angle) < tol) {
74  auto t_R = get_rotation_form_vector(t_omega);
75  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
76  return t_diff_R;
77  }
78 
81  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
82 
83  const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
84  const typename TensorTypeExtractor<T>::Type ss = sin(angle);
85  const typename TensorTypeExtractor<T>::Type a = ss / angle;
86  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
87  const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
88  t_diff_R(i, j, k) = 0;
89  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
90  t_diff_R(i, j, k) +=
91  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
92  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
93 
94  const typename TensorTypeExtractor<T>::Type cc = cos(angle);
95  const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
96  const typename TensorTypeExtractor<T>::Type diff_a =
97  (angle * cc - ss) / angle2;
98 
99  const typename TensorTypeExtractor<T>::Type diff_b =
100  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
101  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
102  t_diff_R(i, j, k) +=
103  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
104 
105  return t_diff_R;
106 }
static Index< 'l', 3 > l
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
double tol
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ get_rotation_form_vector()

template<typename T >
static FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3, 3> EshelbianPlasticity::get_rotation_form_vector ( FTensor::Tensor1< T, 3 > &  t_omega)
static
Examples
EshelbianOperators.cpp.

Definition at line 29 of file EshelbianOperators.cpp.

30  {
35  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
36  t_R(i, j) = t_kd(i, j);
37 
38  const typename TensorTypeExtractor<T>::Type angle =
39  sqrt(t_omega(i) * t_omega(i));
40 
41  constexpr double tol = 1e-18;
42  if (std::abs(angle) < tol) {
43  return t_R;
44  }
45 
47  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
48  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
49  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
50  const typename TensorTypeExtractor<T>::Type b =
51  2. * ss_2 * ss_2 / (angle * angle);
52  t_R(i, j) += a * t_Omega(i, j);
53  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
54 
55  return t_R;
56 }
Kronecker Delta class.
double tol
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ getFTensor3FromMat()

FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3> EshelbianPlasticity::getFTensor3FromMat ( MatrixDouble &  m)
Examples
EshelbianOperators.cpp.

Definition at line 24 of file EshelbianPlasticity.hpp.

24  {
26  &m(0, 0), &m(1, 0), &m(2, 0), &m(3, 0), &m(4, 0), &m(5, 0), &m(6, 0),
27  &m(7, 0), &m(8, 0), &m(9, 0), &m(10, 0), &m(11, 0), &m(12, 0), &m(13, 0),
28  &m(14, 0), &m(15, 0), &m(16, 0), &m(17, 0), &m(18, 0), &m(19, 0),
29  &m(20, 0), &m(21, 0), &m(22, 0), &m(23, 0), &m(24, 0), &m(25, 0),
30  &m(26, 0));
31 }
static Index< 'm', 3 > m

◆ matrix_exponent_derivative()

static MoFEMErrorCode EshelbianPlasticity::matrix_exponent_derivative ( MatrixDouble log_u,
MatrixDouble diff_u,
std::array< MatrixDouble, 6 > &  diff2_u 
)
static
Examples
EshelbianOperators.cpp.

Definition at line 199 of file EshelbianOperators.cpp.

200  {
206 
207  for (auto ll : {0, 1, 2, 3, 4, 5})
208  diff2_u[ll].resize(81, log_u.size2(), false);
209 
210  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
211  t_diff2_u = {
212 
214  &diff2_u[0](0, 0), diff2_u[0].size2()),
216  &diff2_u[1](0, 0), diff2_u[1].size2()),
218  &diff2_u[2](0, 0), diff2_u[2].size2()),
220  &diff2_u[3](0, 0), diff2_u[3].size2()),
222  &diff2_u[4](0, 0), diff2_u[4].size2()),
224  &diff2_u[5](0, 0), diff2_u[5].size2())
225 
226  };
227 
228  constexpr double h = 1e-18;
229  const std::complex<double> ih = h * 1i;
230 
231  int ss = 0;
232  for (auto ii : {0, 1, 2})
233  for (auto jj = 0; jj <= ii; ++jj, ++ss) {
234 
235  auto t_log_u = getFTensor2SymmetricFromMat<3>(log_u);
236  auto t_diff_u =
238  &diff_u(0, 0), diff_u.size2());
239 
241  tc_h(i, j) = 0;
242  tc_h(ii, jj) = ih;
243 
244  for (int gg = 0; gg != log_u.size2(); ++gg) {
245 
248  FTensor::Tensor4<std::complex<double>, 3, 3, 3, 3> tc_diff_u;
249 
250  tc_log_u(i, j) = t_log_u(i, j) + tc_h(i, j);
251  CHKERR tensor_exponent(tc_log_u, tc_u, tc_diff_u);
252 
253  (t_diff2_u[ss])(i, j, k, l) = 0;
254  for (auto kk : {0, 1, 2})
255  for (auto ll : {0, 1, 2})
256  for (auto mm : {0, 1, 2})
257  for (auto nn : {0, 1, 2}) {
258  (t_diff2_u[ss])(kk, ll, mm, nn) =
259  std::imag(tc_diff_u(kk, ll, mm, nn)) / h;
260  }
261 
262  // for (auto kk : {0, 1, 2})
263  // for (auto ll = 0; ll <= 2; ++ll) {
264 
265  // constexpr double eps = 1e-6;
266  // if (std::abs(t_diff_u(kk, ll, ii, jj) -
267  // std::imag(tc_u(kk, ll)) / h) > eps) {
268  // std::ostringstream ss;
269  // ss << "Error " << kk << " " << ll << " : "
270  // << t_diff_u(kk, ll, ii, jj) << " "
271  // << std::imag(tc_u(kk, ll)) / h << " error "
272  // << t_diff_u(kk, ll, ii, jj) - std::imag(tc_u(kk, ll)) / h
273  // << " "
274  // << t_diff_u(kk, ll, ii, jj) / (std::imag(tc_u(kk, ll)) / h)
275  // << endl;
276  // SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s",
277  // ss.str().c_str());
278  // }
279  // }
280 
281  ++t_log_u;
282  ++t_diff_u;
283  ++(t_diff2_u[ss]);
284  }
285  }
286 
288 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static Index< 'l', 3 > l
static MoFEMErrorCode tensor_exponent(FTensor::Tensor2_symmetric< T, 3 > &t_log_s_u, FTensor::Tensor2_symmetric< T, 3 > &t_s_u, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff_s_u)
static Index< 'i', 3 > i
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
static Index< 'k', 3 > k
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

◆ matrix_rotation_derivative()

static MoFEMErrorCode EshelbianPlasticity::matrix_rotation_derivative ( MatrixDouble omega,
std::array< MatrixDouble, 3 > &  diff2_omega 
)
static
Examples
EshelbianOperators.cpp.

Definition at line 291 of file EshelbianOperators.cpp.

292  {
297 
298  for (auto ll : {0, 1, 2})
299  diff2_omega[ll].resize(27, omega.size2(), false);
300 
301  std::array<FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3>, 3>
302  t_diff2_omega = {
303 
304  getFTensor3FromMat(diff2_omega[0]),
305 
306  getFTensor3FromMat(diff2_omega[1]),
307 
308  getFTensor3FromMat(diff2_omega[2])
309 
310  };
311 
312  constexpr double h = 1e-18;
313  const std::complex<double> ih = h * 1i;
314 
315  for (auto ii : {0, 1, 2}) {
316 
317  auto t_omega = getFTensor1FromMat<3>(omega);
318 
320  tc_h(i) = 0;
321  tc_h(ii) = ih;
322 
323  for (int gg = 0; gg != omega.size2(); ++gg) {
324 
326  tc_omega(i) = t_omega(i) + tc_h(i);
327 
328  auto tc_diff_R = get_diff_rotation_form_vector(tc_omega);
329 
330  (t_diff2_omega[ii])(i, j, k) = 0;
331  for (auto kk : {0, 1, 2})
332  for (auto ll : {0, 1, 2})
333  for (auto mm : {0, 1, 2})
334  (t_diff2_omega[ii])(kk, ll, mm) =
335  std::imag(tc_diff_R(kk, ll, mm)) / h;
336 
337  ++t_omega;
338  ++(t_diff2_omega[ii]);
339  }
340  }
342 }
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
static Index< 'i', 3 > i
static Index< 'j', 3 > j
constexpr double omega
static Index< 'k', 3 > k
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

◆ tensor_exponent()

template<typename T >
static MoFEMErrorCode EshelbianPlasticity::tensor_exponent ( FTensor::Tensor2_symmetric< T, 3 > &  t_log_s_u,
FTensor::Tensor2_symmetric< T, 3 > &  t_s_u,
FTensor::Tensor4< T, 3, 3, 3, 3 > &  t_diff_s_u 
)
static
Examples
EshelbianOperators.cpp.

Definition at line 110 of file EshelbianOperators.cpp.

112  {
114 
115  typedef typename TensorTypeExtractor<T>::Type A;
116 
117  auto calc_exponent = [](FTensor::Tensor2<A, 3, 3> &t_log_u,
121 
122  constexpr int N = 800;
123  constexpr double eps = 1e-12;
124 
130 
135 
137  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
138  t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
139 
140  t_p(i, j) = t_log_u(i, j);
141  t_diff_p(i, j, k, l) = t_one4(i, j, k, l);
142 
143  t_u(i, j) = t_p(i, j);
144  for (auto ii : {0, 1, 2})
145  t_u(ii, ii) += 1;
146  t_diff_u(i, j, k, l) = t_diff_p(i, j, k, l);
147 
148  for (int nn = 2; nn <= N; ++nn) {
149 
150  t_diff_tmp(i, j, k, l) =
151  (t_p(i, m) / static_cast<double>(nn)) * t_one4(m, j, k, l) +
152  t_diff_p(i, m, k, l) * (t_log_u(m, j) / static_cast<double>(nn));
153 
154  t_tmp(i, j) = t_p(i, k) * (t_log_u(k, j) / static_cast<double>(nn));
155 
156  t_p(i, j) = t_tmp(i, j);
157  t_diff_p(i, j, k, l) = t_diff_tmp(i, j, k, l);
158 
159  t_u(i, j) += t_p(i, j);
160  t_diff_u(i, j, k, l) += t_diff_p(i, j, k, l);
161 
162  double e = std::abs(sqrt(t_p(i, j) * t_p(i, j)));
163  if (e < eps)
165  }
167  };
168 
171 
173  for (auto ii : {0, 1, 2})
174  for (auto jj : {0, 1, 2})
175  t_log_u(ii, jj) = t_log_s_u(ii, jj);
176 
177  CHKERR calc_exponent(t_log_u, t_u, t_diff_u);
178 
179  for (auto ii : {0, 1, 2})
180  for (auto jj = 0; jj != 3; ++jj) {
181  t_s_u(ii, jj) = (t_u(ii, jj) + t_u(jj, ii)) / 2.;
182  for (auto kk : {0, 1, 2})
183  for (auto ll = 0; ll <= kk; ++ll) {
184 
185  if (ll != kk)
186  t_diff_s_u(ii, jj, kk, ll) =
187 
188  t_diff_u(ii, jj, kk, ll) + t_diff_u(ii, jj, ll, kk);
189 
190  else
191  t_diff_s_u(ii, jj, kk, ll) = t_diff_u(ii, jj, kk, ll);
192  }
193  }
194 
196 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
Kronecker Delta class.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
static Index< 'm', 3 > m
static Index< 'i', 3 > i
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
static Index< 'k', 3 > k
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
static const double eps
const int N
Definition: speed_test.cpp:3

Variable Documentation

◆ IDD_MOFEMEshelbianCrackInterface

const MOFEMuuid EshelbianPlasticity::IDD_MOFEMEshelbianCrackInterface
static
Initial value:
Examples
EshelbianPlasticity.cpp.

Definition at line 17 of file EshelbianPlasticity.hpp.