v0.14.0
Public Member Functions | Private Attributes | List of all members
EshelbianPlasticity::OpCalculateRotationAndSpatialGradient Struct Reference

#include <users_modules/eshelbian_plasticit/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::OpCalculateRotationAndSpatialGradient:
[legend]
Collaboration diagram for EshelbianPlasticity::OpCalculateRotationAndSpatialGradient:
[legend]

Public Member Functions

 OpCalculateRotationAndSpatialGradient (boost::shared_ptr< DataAtIntegrationPts > data_ptr)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 

Private Attributes

boost::shared_ptr< DataAtIntegrationPts > dataAtPts
 data at integration pts More...
 

Detailed Description

Definition at line 557 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpCalculateRotationAndSpatialGradient()

EshelbianPlasticity::OpCalculateRotationAndSpatialGradient::OpCalculateRotationAndSpatialGradient ( boost::shared_ptr< DataAtIntegrationPts >  data_ptr)
inline

Definition at line 559 of file EshelbianPlasticity.hpp.

561  : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}

Member Function Documentation

◆ doWork()

MoFEMErrorCode EshelbianPlasticity::OpCalculateRotationAndSpatialGradient::doWork ( int  side,
EntityType  type,
EntData &  data 
)

Definition at line 231 of file EshelbianOperators.cpp.

233  {
234  MoFEMFunctionBegin;
235 
236  FTensor::Index<'L', size_symm> L;
237  auto t_L = symm_L_tensor();
238 
239  int nb_integration_pts = getGaussPts().size2();
240  i_FTIndex<SPACE_DIM> i;
241  j_FTIndex<SPACE_DIM> j;
242  k_FTIndex<SPACE_DIM> k;
243  l_FTIndex<SPACE_DIM> l;
244  m_FTIndex<SPACE_DIM> m;
245  n_FTIndex<SPACE_DIM> n;
246 
247  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
248  dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
249  dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
250  dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
251  dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
252  dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts, false);
253 
254  dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
255  dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
256  dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
257  dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
258  false);
259  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
260  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
261  dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
262  dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
263  dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
264  dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
265  dataAtPts->nbUniq.resize(nb_integration_pts, false);
266 
267  dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
268  dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
269 
270  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
271  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
272  auto t_h_dlog_u =
273  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
274  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
275  auto t_levi_kirchoff_dP =
276  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
277  auto t_levi_kirchoff_domega =
278  getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
279  auto t_approx_P_adjont_dstretch =
280  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
281  auto t_approx_P_adjont_log_du =
282  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
283  auto t_approx_P_adjont_log_du_dP =
284  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
285  auto t_approx_P_adjont_log_du_domega =
286  getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
287  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
288  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
289  auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
290  auto t_log_u =
291  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
292  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
293  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
294  auto t_diff_u =
295  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
296  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
297  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
298 
299  auto &nbUniq = dataAtPts->nbUniq;
300  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
301  auto t_log_stretch_total =
302  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
303  auto t_log_u2_h1 =
304  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
305  auto t_kd = FTensor::Kronecker_Delta<int>();
306 
307  for (int gg = 0; gg != nb_integration_pts; ++gg) {
308 
309  FTensor::Tensor1<double, 3> eig;
310  FTensor::Tensor2<double, 3, 3> eigen_vec;
311  FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
312  FTensor::Tensor2<double, 3, 3> t_h1;
313  FTensor::Tensor2<double, 3, 3> t_Ru;
314  FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
315 
316  switch (EshelbianCore::gradApproximator) {
317  case NO_H1_CONFIGURATION:
318  t_h1(i, j) = t_kd(i, j);
319  break;
320  case LARGE_ROT:
321  case MODERATE_ROT:
322  t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
323  break;
324  case SMALL_ROT:
325  t_h1(i, j) = t_kd(i, j);
326  break;
327  };
328 
329  auto calculate_rotation = [&]() {
330  auto t0_diff =
331  get_diff_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
332  auto t0 = get_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
333  t_diff_R(i, j, k) = t0_diff(i, j, k);
334  t_R(i, j) = t0(i, j);
335  };
336 
337  auto calculate_stretch = [&]() {
338  if (EshelbianCore::stretchSelector == StretchSelector::LINEAR) {
339  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
340  t_u(i, j) = t_log_u(i, j) + t_kd_sym(i, j);
341  t_diff_u(i, j, k, l) = (t_kd_sym(i, k) ^ t_kd_sym(j, l)) / 4.;
342  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
343  } else {
344  eigen_vec(i, j) = t_log_u(i, j);
345  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
346  // rare case when two eigen values are equal
347  nbUniq[gg] = get_uniq_nb<3>(&eig(0));
348  if (nbUniq[gg] < 3) {
349  sort_eigen_vals<3>(eig, eigen_vec);
350  }
351  t_eigen_vals(i) = eig(i);
352  t_eigen_vecs(i, j) = eigen_vec(i, j);
353  auto t_u_tmp =
354  EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
355  t_u(i, j) = t_u_tmp(i, j);
356  auto t_diff_u_tmp = EigenMatrix::getDiffMat(
357  t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
358  nbUniq[gg]);
359  t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
360  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
361  }
362  };
363 
364  calculate_rotation();
365 
366  if (!EshelbianCore::noStretch) {
367 
368  calculate_stretch();
369 
370  switch (EshelbianCore::gradApproximator) {
371  case NO_H1_CONFIGURATION:
372  case LARGE_ROT:
373 
374  // rotation can be large, moderate or small
375  t_Ru(i, m) = t_R(i, l) * t_u(l, m);
376  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
377  t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
378  t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
379 
380  // conservation of angular momentum at current configuration
381  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
382  t_approx_P_adjont_dstretch(l, m) =
383  t_approx_P_intermidiate(i, m) * t_R(i, l);
384 
385  t_levi_kirchoff(k) =
386  levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
387  t_levi_kirchoff_dP(i, j, k) =
388  (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
389  t_levi_kirchoff_domega(k, n) =
390  levi_civita(l, m, k) *
391  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
392 
393  t_approx_P_adjont_log_du(L) =
394  t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
395  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
396  t_approx_P_adjont_log_du_domega(n, L) =
397  t_Ldiff_u(l, m, L) *
398  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
399 
400  break;
401  case MODERATE_ROT:
402 
403  {
404 
405  // assume small current rotation
406  FTensor::Tensor2<double, 3, 3> t_Omega;
407  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
408  t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
409  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
410  t_h_domega(i, j, k) =
411  FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
412  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
413 
414  // conservation of angular momentum at intermediate configuration
415  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
416  t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
417 
418  t_levi_kirchoff(k) =
419  levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
420  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
421  t_levi_kirchoff_domega(k, n) = 0;
422 
423  t_approx_P_adjont_log_du(L) =
424  t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
425  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
426  t_approx_P_adjont_log_du_domega(n, L) = 0;
427 
428  }
429 
430  break;
431 
432  case SMALL_ROT:
433 
434  {
435 
436  // assume small current rotation
437  FTensor::Tensor2<double, 3, 3> t_Omega;
438  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
439  t_h(i, j) = t_Omega(i, j) + t_u(i, j);
440  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
441  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
442 
443  // conservation of angular momentum at reference condition
444  t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
445  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
446  t_levi_kirchoff_domega(k, l) = 0;
447 
448  t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
449  t_approx_P_adjont_log_du(L) =
450  t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
451  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
452  t_approx_P_adjont_log_du_domega(m, L) = 0;
453 
454  }
455 
456  break;
457  }
458 
459  FTensor::Tensor2<double, 3, 3> t_C_h1;
460  t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
461  eigen_vec(i, j) = t_C_h1(i, j);
462  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
463  switch (EshelbianCore::stretchSelector) {
464  case StretchSelector::LINEAR:
465  break;
466  case StretchSelector::LOG:
467  // for (int ii = 0; ii != 3; ++ii) {
468  // eig(ii) = std::max(
469  // eig(ii), 2 * EshelbianCore::inv_f(EshelbianCore::f(
470  // std::numeric_limits<double>::min_exponent)));
471  // }
472  break;
473  }
474 
475  switch (EshelbianCore::gradApproximator) {
476  case NO_H1_CONFIGURATION:
477  t_log_u2_h1(i, j) = 0;
478  t_log_stretch_total(i, j) = t_log_u(i, j);
479  break;
480  case LARGE_ROT:
481  case MODERATE_ROT: {
482  auto t_log_u2_h1_tmp =
483  EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
484  t_log_u2_h1(i, j) = t_log_u2_h1_tmp(i, j);
485  t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
486  } break;
487  case SMALL_ROT:
488  t_log_u2_h1(i, j) = 0;
489  t_log_stretch_total(i, j) = t_log_u(i, j);
490  break;
491  };
492 
493  } else {
494 
495  auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
496 
497  FTensor::Tensor2<double, 3, 3> t_Omega;
498  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
499  t_u(i, j) = t_kd(i, j) + t_log_u(i, j);
500  t_h(i, j) = t_Omega(i, j) + t_u(i, j);
501  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
502 
503  // conservation of angular momentum at reference condition
504  t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
505  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
506  t_levi_kirchoff_domega(k, l) = 0;
507 
508  t_log_stretch_total(i, j) = t_log_u(i, j);
509  }
510 
511  ++t_h;
512  ++t_h_domega;
513  ++t_h_dlog_u;
514  ++t_levi_kirchoff;
515  ++t_levi_kirchoff_dP;
516  ++t_levi_kirchoff_domega;
517  ++t_approx_P_adjont_dstretch;
518  ++t_approx_P_adjont_log_du;
519  ++t_approx_P_adjont_log_du_dP;
520  ++t_approx_P_adjont_log_du_domega;
521  ++t_approx_P;
522  ++t_R;
523  ++t_diff_R;
524  ++t_log_u;
525  ++t_u;
526  ++t_diff_u;
527  ++t_eigen_vals;
528  ++t_eigen_vecs;
529  ++t_omega;
530  ++t_grad_h1;
531  ++t_log_u2_h1;
532  ++t_log_stretch_total;
533  }
534 
535  MoFEMFunctionReturn(0);
536 }

Member Data Documentation

◆ dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianPlasticity::OpCalculateRotationAndSpatialGradient::dataAtPts
private

data at integration pts

Definition at line 566 of file EshelbianPlasticity.hpp.

The documentation for this struct was generated from the following files:
NOSPACE
@ NOSPACE
Definition: definitions.h:83
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
EshelbianPlasticity::get_rotation_form_vector
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:24
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:44
EshelbianPlasticity::size_symm
constexpr static auto size_symm
Definition: EshelbianAux.hpp:39
FTensor::Tensor1< double, 3 >
EshelbianPlasticity::OpCalculateRotationAndSpatialGradient::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
Definition: EshelbianPlasticity.hpp:566
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity::EshelbianCore::d_f
static boost::function< double(const double)> d_f
Definition: EshelbianPlasticity.hpp:904
EshelbianPlasticity::NO_H1_CONFIGURATION
@ NO_H1_CONFIGURATION
Definition: EshelbianPlasticity.hpp:43
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
FTensor::Tensor2< double, 3, 3 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
EshelbianPlasticity::EshelbianCore::gradApproximator
static enum RotSelector gradApproximator
Definition: EshelbianPlasticity.hpp:897
EshelbianPlasticity::EshelbianCore::stretchSelector
static enum StretchSelector stretchSelector
Definition: EshelbianPlasticity.hpp:898
EshelbianPlasticity::EshelbianCore::rotSelector
static enum RotSelector rotSelector
Definition: EshelbianPlasticity.hpp:896
EshelbianPlasticity::EshelbianCore::f
static boost::function< double(const double)> f
Definition: EshelbianPlasticity.hpp:903
EshelbianPlasticity::EshelbianCore::noStretch
static PetscBool noStretch
Definition: EshelbianPlasticity.hpp:899
EshelbianPlasticity::VolUserDataOperator
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition: EshelbianPlasticity.hpp:51
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
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
EshelbianPlasticity::symm_L_tensor
auto symm_L_tensor()
Definition: EshelbianAux.hpp:52
FTensor::Index
Definition: Index.hpp:23
convert.n
n
Definition: convert.py:82
EshelbianPlasticity::EshelbianCore::inv_f
static boost::function< double(const double)> inv_f
Definition: EshelbianPlasticity.hpp:906
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:44
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
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:60
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:43
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1452
Generated by Doxygen 1.8.17 and hosted at University of Glasgow