v0.14.0
Classes | Functions | Variables
HenckyOps Namespace Reference

Classes

struct  CommonData
 
struct  HenkyIntegrators
 
struct  isEq
 
struct  OpCalculateEigenValsImpl
 
struct  OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpCalculateHenckyPlasticStressImpl
 
struct  OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpCalculateHenckyStressImpl
 
struct  OpCalculateHenckyStressImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpCalculateLogC_dCImpl
 
struct  OpCalculateLogC_dCImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpCalculateLogCImpl
 
struct  OpCalculateLogCImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpCalculatePiolaStressImpl
 
struct  OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >
 
struct  OpHenckyTangentImpl
 
struct  OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >
 

Functions

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)
 
template<int DIM>
MoFEMErrorCode addMatBlockOps (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
 
template<int DIM, IntegrationType I, typename DomainEleOp >
auto commonDataFactory (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
 

Variables

constexpr double eps = std::numeric_limits<float>::epsilon()
 
auto f = [](double v) { return 0.5 * std::log(v); }
 
auto d_f = [](double v) { return 0.5 / v; }
 
auto dd_f = [](double v) { return -0.5 / (v * v); }
 

Function Documentation

◆ addMatBlockOps()

template<int DIM>
MoFEMErrorCode HenckyOps::addMatBlockOps ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  block_name,
boost::shared_ptr< MatrixDouble >  mat_D_Ptr,
Sev  sev,
double  scale = 1 
)

[Calculate elasticity tensor]

[Calculate elasticity tensor]

Definition at line 597 of file HenckyOps.hpp.

601  {
602  MoFEMFunctionBegin;
603 
604  struct OpMatBlocks : public DomainEleOp {
605  OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
606  double shear_modulus_G, MoFEM::Interface &m_field, Sev sev,
607  std::vector<const CubitMeshSets *> meshset_vec_ptr,
608  double scale)
609  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
610  bulkModulusKDefault(bulk_modulus_K),
611  shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
612  CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
613  "Can not get data from block");
614  }
615 
616  MoFEMErrorCode doWork(int side, EntityType type,
617  EntitiesFieldData::EntData &data) {
618  MoFEMFunctionBegin;
619 
620  for (auto &b : blockData) {
621 
622  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
623  CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
624  b.shearModulusG * scaleYoungModulus);
625  MoFEMFunctionReturnHot(0);
626  }
627  }
628 
629  CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
630  shearModulusGDefault * scaleYoungModulus);
631  MoFEMFunctionReturn(0);
632  }
633 
634  private:
635  boost::shared_ptr<MatrixDouble> matDPtr;
636  const double scaleYoungModulus;
637 
638  struct BlockData {
639  double bulkModulusK;
640  double shearModulusG;
641  Range blockEnts;
642  };
643 
644  double bulkModulusKDefault;
645  double shearModulusGDefault;
646  std::vector<BlockData> blockData;
647 
648  MoFEMErrorCode
649  extractBlockData(MoFEM::Interface &m_field,
650  std::vector<const CubitMeshSets *> meshset_vec_ptr,
651  Sev sev) {
652  MoFEMFunctionBegin;
653 
654  for (auto m : meshset_vec_ptr) {
655  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
656  std::vector<double> block_data;
657  CHKERR m->getAttributes(block_data);
658  if (block_data.size() != 2) {
659  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
660  "Expected that block has two attribute");
661  }
662  auto get_block_ents = [&]() {
663  Range ents;
664  CHKERR
665  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
666  return ents;
667  };
668 
669  double young_modulus = block_data[0];
670  double poisson_ratio = block_data[1];
671  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
672  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
673 
674  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
675  << "E = " << young_modulus << " nu = " << poisson_ratio;
676 
677  blockData.push_back(
678  {bulk_modulus_K, shear_modulus_G, get_block_ents()});
679  }
680  MOFEM_LOG_CHANNEL("WORLD");
681  MoFEMFunctionReturn(0);
682  }
683 
684  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
685  double bulk_modulus_K, double shear_modulus_G) {
686  MoFEMFunctionBegin;
687  //! [Calculate elasticity tensor]
688  auto set_material_stiffness = [&]() {
689  FTensor::Index<'i', DIM> i;
690  FTensor::Index<'j', DIM> j;
691  FTensor::Index<'k', DIM> k;
692  FTensor::Index<'l', DIM> l;
693  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
694  double A = (DIM == 2)
695  ? 2 * shear_modulus_G /
696  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
697  : 1;
698  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
699  t_D(i, j, k, l) =
700  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
701  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
702  t_kd(k, l);
703  };
704  //! [Calculate elasticity tensor]
705  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
706  mat_D_ptr->resize(size_symm * size_symm, 1);
707  set_material_stiffness();
708  MoFEMFunctionReturn(0);
709  }
710  };
711 
712  double E = young_modulus;
713  double nu = poisson_ratio;
714 
715  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
716  CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
717  PETSC_NULL);
718  CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
719  PETSC_NULL);
720  ierr = PetscOptionsEnd();
721 
722  double bulk_modulus_K = E / (3 * (1 - 2 * nu));
723  double shear_modulus_G = E / (2 * (1 + nu));
724  pip.push_back(new OpMatBlocks(
725  mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
726 
727  // Get blockset using regular expression
728  m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
729 
730  (boost::format("%s(.*)") % block_name).str()
731 
732  )),
733  scale
734 
735  ));
736 
737  MoFEMFunctionReturn(0);
738 }

◆ commonDataFactory()

template<int DIM, IntegrationType I, typename DomainEleOp >
auto HenckyOps::commonDataFactory ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  field_name,
std::string  block_name,
Sev  sev,
double  scale = 1 
)

Definition at line 741 of file HenckyOps.hpp.

744  {
745 
746  auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
747  common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
748  common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
749 
750  CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
751  common_ptr->matDPtr, sev, scale),
752  "addMatBlockOps");
753 
754  using H = HenkyIntegrators<DomainEleOp>;
755 
756  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
757  field_name, common_ptr->matGradPtr));
758  pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
759  field_name, common_ptr));
760  pip.push_back(
761  new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
762  pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
763  field_name, common_ptr));
764  pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I>(
765  field_name, common_ptr));
766  pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
767  field_name, common_ptr));
768 
769  return common_ptr;
770 }

◆ get_uniq_nb()

template<int DIM>
auto HenckyOps::get_uniq_nb ( double *  ptr)
inline

Definition at line 32 of file HenckyOps.hpp.

32  {
33  std::array<double, DIM> tmp;
34  std::copy(ptr, &ptr[DIM], tmp.begin());
35  std::sort(tmp.begin(), tmp.end());
36  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
37  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
38 };

◆ is_eq()

auto HenckyOps::is_eq ( const double &  a,
const double &  b 
)
inline

Definition at line 28 of file HenckyOps.hpp.

28  {
29  return isEq::check(a, b);
30 };

◆ opFactoryDomainLhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode HenckyOps::opFactoryDomainLhs ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  field_name,
boost::shared_ptr< HenckyOps::CommonData >  common_ptr,
Sev  sev 
)

Definition at line 806 of file HenckyOps.hpp.

810  {
811  MoFEMFunctionBegin;
812 
813  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
814  A>::template BiLinearForm<I>;
815  using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
816 
817  using H = HenkyIntegrators<DomainEleOp>;
818  pip.push_back(
819  new typename H::template OpHenckyTangent<DIM, I>(field_name, common_ptr));
820  pip.push_back(
821  new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
822 
823  MoFEMFunctionReturn(0);
824 }

◆ opFactoryDomainLhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode HenckyOps::opFactoryDomainLhs ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  field_name,
std::string  block_name,
Sev  sev,
double  scale = 1 
)

Definition at line 827 of file HenckyOps.hpp.

830  {
831  MoFEMFunctionBegin;
832 
833  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
834  m_field, pip, field_name, block_name, sev, scale);
835  CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
836  common_ptr, sev);
837 
838  MoFEMFunctionReturn(0);
839 }

◆ opFactoryDomainRhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode HenckyOps::opFactoryDomainRhs ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  field_name,
boost::shared_ptr< HenckyOps::CommonData >  common_ptr,
Sev  sev 
)

Definition at line 773 of file HenckyOps.hpp.

777  {
778  MoFEMFunctionBegin;
779 
780  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
781  A>::template LinearForm<I>;
782  using OpInternalForcePiola =
783  typename B::template OpGradTimesTensor<1, DIM, DIM>;
784  pip.push_back(
785  new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
786 
787  MoFEMFunctionReturn(0);
788 }

◆ opFactoryDomainRhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode HenckyOps::opFactoryDomainRhs ( MoFEM::Interface &  m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pip,
std::string  field_name,
std::string  block_name,
Sev  sev,
double  scale = 1 
)

Definition at line 791 of file HenckyOps.hpp.

794  {
795  MoFEMFunctionBegin;
796 
797  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
798  m_field, pip, field_name, block_name, sev, scale);
799  CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
800  common_ptr, sev);
801 
802  MoFEMFunctionReturn(0);
803 }

◆ sort_eigen_vals()

template<int DIM>
auto HenckyOps::sort_eigen_vals ( FTensor::Tensor1< double, DIM > &  eig,
FTensor::Tensor2< double, DIM, DIM > &  eigen_vec 
)
inline

Definition at line 41 of file HenckyOps.hpp.

42  {
43 
44  isEq::absMax =
45  std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
46 
47  int i = 0, j = 1, k = 2;
48 
49  if (is_eq(eig(0), eig(1))) {
50  i = 0;
51  j = 2;
52  k = 1;
53  } else if (is_eq(eig(0), eig(2))) {
54  i = 0;
55  j = 1;
56  k = 2;
57  } else if (is_eq(eig(1), eig(2))) {
58  i = 1;
59  j = 0;
60  k = 2;
61  }
62 
63  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
64  eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
65 
66  eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
67 
68  eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
69 
70  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
71 
72  {
73  FTensor::Index<'i', DIM> i;
74  FTensor::Index<'j', DIM> j;
75  eig(i) = eig_c(i);
76  eigen_vec(i, j) = eigen_vec_c(i, j);
77  }
78 };

Variable Documentation

◆ d_f

auto HenckyOps::d_f = [](double v) { return 0.5 / v; }
Examples
EshelbianOperators.cpp, EshelbianPlasticity.cpp, and matrix_function.cpp.

Definition at line 16 of file HenckyOps.hpp.

◆ dd_f

auto HenckyOps::dd_f = [](double v) { return -0.5 / (v * v); }
Examples
EshelbianOperators.cpp, EshelbianPlasticity.cpp, and matrix_function.cpp.

Definition at line 17 of file HenckyOps.hpp.

◆ eps

constexpr double HenckyOps::eps = std::numeric_limits<float>::epsilon()
constexpr

Definition at line 13 of file HenckyOps.hpp.

◆ f

auto HenckyOps::f = [](double v) { return 0.5 * std::log(v); }
Examples
ADOLCPlasticityMaterialModels.hpp, bernstein_bezier_generate_base.cpp, dm_create_subdm.cpp, EshelbianOperators.cpp, EshelbianPlasticity.cpp, free_surface.cpp, matrix_function.cpp, NavierStokesElement.cpp, plastic.cpp, PlasticOpsGeneric.hpp, PoissonOperators.hpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, reaction_diffusion.cpp, scalar_check_approximation.cpp, shallow_wave.cpp, simple_elasticity.cpp, simple_interface.cpp, tensor_divergence_operator.cpp, testing_jacobian_of_hook_element.cpp, testing_jacobian_of_hook_scaled_with_density_element.cpp, and UnsaturatedFlow.hpp.

Definition at line 15 of file HenckyOps.hpp.

NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
FTensor::Tensor1< double, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
FTensor::Tensor2< double, 3, 3 >
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
HenckyOps::is_eq
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:28
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
BlockData
Definition: gel_analysis.cpp:74
a
constexpr double a
Definition: approx_sphere.cpp:30
convert.type
type
Definition: convert.py:64
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
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
BiLinearForm
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', DIM >
Range
DomainEleOp
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
H
double H
Hardening.
Definition: plastic.cpp:175
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
Generated by Doxygen 1.8.17 and hosted at University of Glasgow