v0.9.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  OpPostProcDataStructure
 
struct  OpRotationBc
 
struct  OpRotationBc_dx
 
struct  OpSpatialConsistency_dBubble_domega
 
struct  OpSpatialConsistency_dBubble_dx
 
struct  OpSpatialConsistency_dP_domega
 
struct  OpSpatialConsistency_dP_dx
 
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  OpSpatialRotation_domega_du
 
struct  OpSpatialRotation_domega_dx
 
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_dot_deviator (FTensor::Tensor2< T, 3, 3 > &t)
 
template<typename T >
static FTensor::Tensor4< typename TensorTypeExtractor< T >::Type, 3, 3, 3, 3 > calculate_diff_dot_deviator (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 530 of file EshelbianPlasticity.hpp.

◆ BcRotVec

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

Definition at line 539 of file EshelbianPlasticity.hpp.

◆ EntData

using EshelbianPlasticity::EntData = typedef DataForcesAndSourcesCore::EntData

Definition at line 34 of file EshelbianPlasticity.hpp.

◆ FaceUserDataOperator

using EshelbianPlasticity::FaceUserDataOperator = typedef FaceElementForcesAndSourcesCore::UserDataOperator

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 550 of file EshelbianPlasticity.hpp.

◆ TractionFreeBc

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

Definition at line 541 of file EshelbianPlasticity.hpp.

◆ UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

◆ VectorPtr

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

Definition at line 21 of file EshelbianPlasticity.hpp.

◆ VolUserDataOperator

using EshelbianPlasticity::VolUserDataOperator = typedef VolumeElementForcesAndSourcesCore::UserDataOperator

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_dot_deviator()

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

Definition at line 385 of file EshelbianOperators.cpp.

386  {
387 
393 
394  auto get_dot_fun = [i, j, k, m, n](auto &t, auto &t_diff) {
396  t_diff_c;
397  t_diff_c(i, j, m, n) = t_diff(i, j, m, n);
398  return t_diff_c;
399  };
400 
401  auto get_trace = [i, j](auto &t_diff) {
403  t_diff_tr(i, j) = 0;
404  for (auto ii = 0; ii != 3; ++ii)
405  for (auto mm = 0; mm != 3; ++mm)
406  for (auto nn = 0; nn != 3; ++nn)
407  t_diff_tr(mm, nn) += t_diff(ii, ii, mm, nn) / 3.;
408  return t_diff_tr;
409  };
410 
411  auto get_deviator = [i, j, m, n, get_trace](auto &&t_diff) {
412  auto t_diff_tr = get_trace(t_diff);
414  t_diff_dev;
415  t_diff_dev(i, j, m, n) = t_diff(i, j, m, n);
416  for (auto mm = 0; mm != 3; ++mm)
417  for (auto nn = 0; nn != 3; ++nn)
418  for (auto ii = 0; ii != 3; ++ii)
419  t_diff_dev(ii, ii, mm, nn) -= t_diff_tr(mm, nn);
420  return t_diff_dev;
421  };
422 
423  return get_deviator(get_dot_fun(t, t_diff));
424 }

◆ calculate_dot_deviator()

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

Definition at line 352 of file EshelbianOperators.cpp.

352  {
353 
357 
358  auto get_dot_fun = [i, j, k](auto &t) {
360  t_c(i, j) = t(i, j);
361  return t_c;
362  };
363 
364  auto get_trace = [](auto &t) {
365  typename TensorTypeExtractor<T>::Type tr = 0;
366  for (auto ii = 0; ii != 3; ++ii)
367  tr += t(ii, ii) / 3.;
368  return tr;
369  };
370 
371  auto get_deviator = [i, j, get_trace](auto &&t) {
372  const auto tr = get_trace(t);
374  t_dev(i, j) = t(i, j);
375  for (auto ii = 0; ii != 3; ++ii)
376  t_dev(ii, ii) -= tr;
377  return t_dev;
378  };
379 
380  return get_deviator(get_dot_fun(t));
381 }

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

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 
31  Index<'i', 3> i;
32  Index<'j', 3> j;
34  Index<'m', 3> m;
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:505
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:512
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
#define CHKERR
Inline error check.
Definition: definitions.h:600
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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 62 of file EshelbianOperators.cpp.

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

◆ 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  t_R(i, j) = 0.;
36  t_R(0, 0) = 1.;
37  t_R(1, 1) = 1.;
38  t_R(2, 2) = 1.;
39 
40  const typename TensorTypeExtractor<T>::Type angle =
41  sqrt(t_omega(i) * t_omega(i));
42 
43  constexpr double tol = 1e-18;
44  if (std::abs(angle) < tol) {
45  return t_R;
46  }
47 
49  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
50  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
51  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
52  const typename TensorTypeExtractor<T>::Type b =
53  2. * ss_2 * ss_2 / (angle * angle);
54  t_R(i, j) += a * t_Omega(i, j);
55  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
56 
57  return t_R;
58 }
double tol

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

◆ 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 205 of file EshelbianOperators.cpp.

206  {
212 
213  for (auto ll : {0, 1, 2, 3, 4, 5})
214  diff2_u[ll].resize(81, log_u.size2(), false);
215 
216  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
217  t_diff2_u = {
218 
220  &diff2_u[0](0, 0), diff2_u[0].size2()),
222  &diff2_u[1](0, 0), diff2_u[1].size2()),
224  &diff2_u[2](0, 0), diff2_u[2].size2()),
226  &diff2_u[3](0, 0), diff2_u[3].size2()),
228  &diff2_u[4](0, 0), diff2_u[4].size2()),
230  &diff2_u[5](0, 0), diff2_u[5].size2())
231 
232  };
233 
234  constexpr double h = 1e-18;
235  const std::complex<double> ih = h * 1i;
236 
237  int ss = 0;
238  for (auto ii : {0, 1, 2})
239  for (auto jj = 0; jj <= ii; ++jj, ++ss) {
240 
241  auto t_log_u = getFTensor2SymmetricFromMat<3>(log_u);
242  auto t_diff_u =
244  &diff_u(0, 0), diff_u.size2());
245 
247  tc_h(i, j) = 0;
248  tc_h(ii, jj) = ih;
249 
250  for (int gg = 0; gg != log_u.size2(); ++gg) {
251 
254  FTensor::Tensor4<std::complex<double>, 3, 3, 3, 3> tc_diff_u;
255 
256  tc_log_u(i, j) = t_log_u(i, j) + tc_h(i, j);
257  CHKERR tensor_exponent(tc_log_u, tc_u, tc_diff_u);
258 
259  (t_diff2_u[ss])(i, j, k, l) = 0;
260  for (auto kk : {0, 1, 2})
261  for (auto ll : {0, 1, 2})
262  for (auto mm : {0, 1, 2})
263  for (auto nn : {0, 1, 2}) {
264  (t_diff2_u[ss])(kk, ll, mm, nn) =
265  std::imag(tc_diff_u(kk, ll, mm, nn)) / h;
266  }
267 
268  // for (auto kk : {0, 1, 2})
269  // for (auto ll = 0; ll <= 2; ++ll) {
270 
271  // constexpr double eps = 1e-6;
272  // if (std::abs(t_diff_u(kk, ll, ii, jj) -
273  // std::imag(tc_u(kk, ll)) / h) > eps) {
274  // std::ostringstream ss;
275  // ss << "Error " << kk << " " << ll << " : "
276  // << t_diff_u(kk, ll, ii, jj) << " "
277  // << std::imag(tc_u(kk, ll)) / h << " error "
278  // << t_diff_u(kk, ll, ii, jj) - std::imag(tc_u(kk, ll)) / h
279  // << " "
280  // << t_diff_u(kk, ll, ii, jj) / (std::imag(tc_u(kk, ll)) / h)
281  // << endl;
282  // SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s",
283  // ss.str().c_str());
284  // }
285  // }
286 
287  ++t_log_u;
288  ++t_diff_u;
289  ++(t_diff2_u[ss]);
290  }
291  }
292 
294 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
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)
#define CHKERR
Inline error check.
Definition: definitions.h:600
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411

◆ matrix_rotation_derivative()

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

Definition at line 297 of file EshelbianOperators.cpp.

298  {
303 
304  for (auto ll : {0, 1, 2})
305  diff2_omega[ll].resize(27, omega.size2(), false);
306 
307  std::array<FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3>, 3>
308  t_diff2_omega = {
309 
310  getFTensor3FromMat(diff2_omega[0]),
311 
312  getFTensor3FromMat(diff2_omega[1]),
313 
314  getFTensor3FromMat(diff2_omega[2])
315 
316  };
317 
318  constexpr double h = 1e-18;
319  const std::complex<double> ih = h * 1i;
320 
321  for (auto ii : {0, 1, 2}) {
322 
323  auto t_omega = getFTensor1FromMat<3>(omega);
324 
326  tc_h(i) = 0;
327  tc_h(ii) = ih;
328 
329  for (int gg = 0; gg != omega.size2(); ++gg) {
330 
332  tc_omega(i) = t_omega(i) + tc_h(i);
333 
334  auto tc_diff_R = get_diff_rotation_form_vector(tc_omega);
335 
336  (t_diff2_omega[ii])(i, j, k) = 0;
337  for (auto kk : {0, 1, 2})
338  for (auto ll : {0, 1, 2})
339  for (auto mm : {0, 1, 2})
340  (t_diff2_omega[ii])(kk, ll, mm) =
341  std::imag(tc_diff_R(kk, ll, mm)) / h;
342 
343  ++t_omega;
344  ++(t_diff2_omega[ii]);
345  }
346  }
348 }
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:481
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411

◆ 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 112 of file EshelbianOperators.cpp.

114  {
116 
117  typedef typename TensorTypeExtractor<T>::Type A;
118 
119  auto calc_exponent = [](FTensor::Tensor2<A, 3, 3> &t_log_u,
123 
124  constexpr int N = 800;
125  constexpr double eps = 1e-12;
126 
132 
137 
140 
141  t_one2(i, j) = 0;
142  for (auto ii : {0, 1, 2})
143  t_one2(ii, ii) = 1;
144  t_one4(i, j, k, l) = t_one2(j, l) * t_one2(i, k);
145 
146  t_p(i, j) = t_log_u(i, j);
147  t_diff_p(i, j, k, l) = t_one4(i, j, k, l);
148 
149  t_u(i, j) = t_p(i, j);
150  for (auto ii : {0, 1, 2})
151  t_u(ii, ii) += 1;
152  t_diff_u(i, j, k, l) = t_diff_p(i, j, k, l);
153 
154  for (int nn = 2; nn <= N; ++nn) {
155 
156  t_diff_tmp(i, j, k, l) =
157  (t_p(i, m) / static_cast<double>(nn)) * t_one4(m, j, k, l) +
158  t_diff_p(i, m, k, l) * (t_log_u(m, j) / static_cast<double>(nn));
159 
160  t_tmp(i, j) = t_p(i, k) * (t_log_u(k, j) / static_cast<double>(nn));
161 
162  t_p(i, j) = t_tmp(i, j);
163  t_diff_p(i, j, k, l) = t_diff_tmp(i, j, k, l);
164 
165  t_u(i, j) += t_p(i, j);
166  t_diff_u(i, j, k, l) += t_diff_p(i, j, k, l);
167 
168  double e = std::abs(sqrt(t_p(i, j) * t_p(i, j)));
169  if (e < eps)
171  }
173  };
174 
177 
179  for (auto ii : {0, 1, 2})
180  for (auto jj : {0, 1, 2})
181  t_log_u(ii, jj) = t_log_s_u(ii, jj);
182 
183  CHKERR calc_exponent(t_log_u, t_u, t_diff_u);
184 
185  for (auto ii : {0, 1, 2})
186  for (auto jj = 0; jj != 3; ++jj) {
187  t_s_u(ii, jj) = (t_u(ii, jj) + t_u(jj, ii)) / 2.;
188  for (auto kk : {0, 1, 2})
189  for (auto ll = 0; ll <= kk; ++ll) {
190 
191  if (ll != kk)
192  t_diff_s_u(ii, jj, kk, ll) =
193 
194  t_diff_u(ii, jj, kk, ll) + t_diff_u(ii, jj, ll, kk);
195 
196  else
197  t_diff_s_u(ii, jj, kk, ll) = t_diff_u(ii, jj, kk, ll);
198  }
199  }
200 
202 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:505
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:512
#define CHKERR
Inline error check.
Definition: definitions.h:600
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411
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.