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

Classes

struct  BcDisp
 
struct  BcLoadScale
 
struct  BcRot
 
struct  CGGUserPolynomialBase
 
struct  DataAtIntegrationPts
 
struct  EshelbianCore
 
struct  FaceRule
 
struct  FeTractionFreeBc
 
struct  HMHStVenantKirchhoff
 
struct  OpAssembleBasic
 
struct  OpAssembleFace
 
struct  OpAssembleVolume
 
struct  OpCalculateRotationAndSpatialGradient
 
struct  OpDispBc
 
struct  OpHMHHStVenantKirchhoff
 
struct  OpJacobian
 
struct  OpPostProcDataStructure
 
struct  OpRotationBc
 
struct  OpSpatialConsistency_dBubble_domega
 
struct  OpSpatialConsistency_dP_domega
 
struct  OpSpatialConsistencyBubble
 
struct  OpSpatialConsistencyDivTerm
 
struct  OpSpatialConsistencyP
 
struct  OpSpatialEquilibrium
 
struct  OpSpatialEquilibrium_dw_dP
 
struct  OpSpatialPhysical
 
struct  OpSpatialPhysical_du_dBubble
 
struct  OpSpatialPhysical_du_domega
 
struct  OpSpatialPhysical_du_dP
 
struct  OpSpatialPhysical_du_du
 
struct  OpSpatialRotation
 
struct  OpSpatialRotation_domega_dBubble
 
struct  OpSpatialRotation_domega_domega
 
struct  OpSpatialRotation_domega_dP
 
struct  PhysicalEquations
 
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
 

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 >
FTensor::Tensor2< double, 3, 3 > getRotationFormVector (FTensor::Tensor1< T, 3 > &t_omega)
 
template<typename T >
FTensor::Tensor3< double, 3, 3, 3 > getDiffRotationFormVector (FTensor::Tensor1< T, 3 > &t_omega)
 

Variables

static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
 

Typedef Documentation

◆ BcDispVec

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

Definition at line 337 of file EshelbianPlasticity.hpp.

◆ BcRotVec

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

Definition at line 345 of file EshelbianPlasticity.hpp.

◆ EntData

using EshelbianPlasticity::EntData = typedef DataForcesAndSourcesCore::EntData

◆ FaceUserDataOperator

using EshelbianPlasticity::FaceUserDataOperator = typedef FaceElementForcesAndSourcesCore::UserDataOperator

Definition at line 107 of file EshelbianPlasticity.hpp.

◆ MatrixPtr

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

Definition at line 16 of file EshelbianPlasticity.hpp.

◆ TractionFreeBc

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

Definition at line 347 of file EshelbianPlasticity.hpp.

◆ UserDataOperator

using EshelbianPlasticity::UserDataOperator = typedef ForcesAndSourcesCore::UserDataOperator

◆ VectorPtr

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

Definition at line 17 of file EshelbianPlasticity.hpp.

◆ VolUserDataOperator

using EshelbianPlasticity::VolUserDataOperator = typedef VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 106 of file EshelbianPlasticity.hpp.

Enumeration Type Documentation

◆ EshelbianInterfaces

Enumerator
CORE_ESHELBIAN_INTERFACE 

Definition at line 11 of file EshelbianPlasticity.hpp.

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 [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

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;
33  Index<'k', 3> k;
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  const double l2 = l_i * l_j * l_k;
153  Tensor1<double, 3> t_diff_l2;
154  t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
155  t_diff_k(i) * l_i * l_j;
156  Tensor2<double, 3, 3> t_diff2_l2;
157  t_diff2_l2(i, j) =
158  t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
159  t_diff_i(i) * l_j * t_diff_k(j) +
160 
161  t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
162  t_diff_j(i) * l_i * t_diff_k(j) +
163 
164  t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
165  t_diff_k(i) * l_i * t_diff_j(j);
166 
167  for (int dd = 0; dd != 3; ++dd) {
168 
169  Tensor2<double, 3, 3> t_axial_diff;
170  t_axial_diff(i, j) = 0;
171  for (int mm = 0; mm != 3; ++mm)
172  t_axial_diff(dd, mm) = t_diff_l2(mm);
173 
174  Tensor3<double, 3, 3, 3> t_A_diff;
175  t_A_diff(i, j, k) = levi_civita(i, j, m) * t_axial_diff(m, k);
176  Tensor2<double, 3, 3> t_curl_A;
177  t_curl_A(i, j) = levi_civita(j, m, f) * t_A_diff(i, f, m);
178  Tensor3<double, 3, 3, 3> t_curl_A_bK_diff;
179  t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
180 
181  Tensor3<double, 3, 3, 3> t_axial_diff2;
182  t_axial_diff2(i, j, k) = 0;
183  for (int mm = 0; mm != 3; ++mm)
184  for (int nn = 0; nn != 3; ++nn)
185  t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
186  Tensor4<double, 3, 3, 3, 3> t_A_diff2;
187  t_A_diff2(i, j, k, f) =
188  levi_civita(i, j, m) * t_axial_diff2(m, k, f);
189  Tensor3<double, 3, 3, 3> t_curl_A_diff2;
190  t_curl_A_diff2(i, j, k) =
191  levi_civita(j, m, f) * t_A_diff2(i, f, m, k);
192  Tensor3<double, 3, 3, 3> t_curl_A_diff2_bK;
193  t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
194 
195  t_phi(i, j) = levi_civita(j, m, f) * (t_curl_A_bK_diff(i, f, m) +
196  t_curl_A_diff2_bK(i, f, m));
197 
198  ++t_phi;
199  ++zz;
200  }
201  }
202  }
203  if(zz != NBVOLUMETET_CCG_BUBBLE(p))
204  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
205  "Wrong number of base functions %d != %d", zz,
207 
208  }
209 
211 }
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
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 MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Common.hpp:212
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Common.hpp:211
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
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

◆ getDiffRotationFormVector()

template<typename T >
FTensor::Tensor3<double, 3, 3, 3> EshelbianPlasticity::getDiffRotationFormVector ( FTensor::Tensor1< T, 3 > &  t_omega)

Definition at line 61 of file EshelbianPlasticity.hpp.

61  {
62 
67 
68  double angle = sqrt(t_omega(i) * t_omega(i));
69  const double tol = 1e-12;
70  if (fabs(angle) < tol) {
72  auto t_R = getRotationFormVector(t_omega);
73  t_diff_R(i, j, k) = levi_civita(i, l, k) * t_R(l, j);
74  return t_diff_R;
75  }
76 
79  t_Omega(i, j) = levi_civita(i, j, k) * t_omega(k);
80 
81  const double angle2 = angle * angle;
82  const double ss = sin(angle);
83  const double a = ss / angle;
84  const double ss_2 = sin(angle / 2.);
85  const double b = 2 * ss_2 * ss_2 / angle2;
86  t_diff_R(i, j, k) = 0;
87  t_diff_R(i, j, k) = a * levi_civita(i, j, k);
88  t_diff_R(i, j, k) += b * (t_Omega(i, l) * levi_civita(l, j, k) +
89  levi_civita(i, l, k) * t_Omega(l, j));
90 
91  const double cc = cos(angle);
92  const double cc_2 = cos(angle / 2.);
93  const double diff_a = (angle * cc - ss) / angle2;
94 
95  const double diff_b =
96  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
97  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
98  t_diff_R(i, j, k) +=
99  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
100 
101  return t_diff_R;
102 }
double tol
FTensor::Tensor2< double, 3, 3 > getRotationFormVector(FTensor::Tensor1< T, 3 > &t_omega)
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

◆ getFTensor3FromMat()

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

Definition at line 20 of file EshelbianPlasticity.hpp.

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

◆ getRotationFormVector()

template<typename T >
FTensor::Tensor2<double, 3, 3> EshelbianPlasticity::getRotationFormVector ( FTensor::Tensor1< T, 3 > &  t_omega)

Definition at line 31 of file EshelbianPlasticity.hpp.

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

Variable Documentation

◆ IDD_MOFEMEshelbianCrackInterface

const MOFEMuuid EshelbianPlasticity::IDD_MOFEMEshelbianCrackInterface
static
Initial value:

Definition at line 13 of file EshelbianPlasticity.hpp.