v0.14.0
Loading...
Searching...
No Matches
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 field_name, 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  field_name,
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 {
603
604 struct OpMatBlocks : public DomainEleOp {
605 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
606 double bulk_modulus_K, double shear_modulus_G,
607 MoFEM::Interface &m_field, Sev sev,
608 std::vector<const CubitMeshSets *> meshset_vec_ptr,
609 double scale)
610 : DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
611 bulkModulusKDefault(bulk_modulus_K),
612 shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
613 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
614 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
615 "Can not get data from block");
616 }
617
618 MoFEMErrorCode doWork(int side, EntityType type,
619 EntitiesFieldData::EntData &data) {
621
622 for (auto &b : blockData) {
623
624 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
625 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
626 b.shearModulusG * scaleYoungModulus);
628 }
629 }
630
631 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
632 shearModulusGDefault * scaleYoungModulus);
634 }
635
636 private:
637 boost::shared_ptr<MatrixDouble> matDPtr;
638 const double scaleYoungModulus;
639
640 struct BlockData {
641 double bulkModulusK;
642 double shearModulusG;
643 Range blockEnts;
644 };
645
646 double bulkModulusKDefault;
647 double shearModulusGDefault;
648 std::vector<BlockData> blockData;
649
651 extractBlockData(MoFEM::Interface &m_field,
652 std::vector<const CubitMeshSets *> meshset_vec_ptr,
653 Sev sev) {
655
656 for (auto m : meshset_vec_ptr) {
657 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
658 std::vector<double> block_data;
659 CHKERR m->getAttributes(block_data);
660 if (block_data.size() != 2) {
661 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
662 "Expected that block has two attribute");
663 }
664 auto get_block_ents = [&]() {
665 Range ents;
666 CHKERR
667 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
668 return ents;
669 };
670
671 double young_modulus = block_data[0];
672 double poisson_ratio = block_data[1];
673 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
674 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
675
676 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
677 << "E = " << young_modulus << " nu = " << poisson_ratio;
678
679 blockData.push_back(
680 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
681 }
682 MOFEM_LOG_CHANNEL("WORLD");
684 }
685
686 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
687 double bulk_modulus_K, double shear_modulus_G) {
689 //! [Calculate elasticity tensor]
690 auto set_material_stiffness = [&]() {
696 double A = (DIM == 2)
697 ? 2 * shear_modulus_G /
698 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
699 : 1;
700 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
701 t_D(i, j, k, l) =
702 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
703 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
704 t_kd(k, l);
705 };
706 //! [Calculate elasticity tensor]
707 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
708 mat_D_ptr->resize(size_symm * size_symm, 1);
709 set_material_stiffness();
711 }
712 };
713
714 double E = young_modulus;
715 double nu = poisson_ratio;
716
717 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
718 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
719 PETSC_NULL);
720 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
721 PETSC_NULL);
722 ierr = PetscOptionsEnd();
723
724 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
725 double shear_modulus_G = E / (2 * (1 + nu));
726 pip.push_back(new OpMatBlocks(
727 field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
728
729 // Get blockset using regular expression
730 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
731
732 (boost::format("%s(.*)") % block_name).str()
733
734 )),
735 scale
736
737 ));
738
740}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
static PetscErrorCode ierr
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
double bulk_modulus_K
double shear_modulus_G
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
constexpr AssemblyType A
double young_modulus
Young modulus.
Definition: plastic.cpp:172
double scale
Definition: plastic.cpp:170
constexpr auto size_symm
Definition: plastic.cpp:42
constexpr auto field_name
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ 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 743 of file HenckyOps.hpp.

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

◆ 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 808 of file HenckyOps.hpp.

812 {
814
815 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
816 A>::template BiLinearForm<I>;
817 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
818
820 pip.push_back(
821 new typename H::template OpHenckyTangent<DIM, I>(field_name, common_ptr));
822 pip.push_back(
823 new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
824
826}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64

◆ 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 829 of file HenckyOps.hpp.

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

◆ 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 775 of file HenckyOps.hpp.

779 {
781
782 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
783 A>::template LinearForm<I>;
785 typename B::template OpGradTimesTensor<1, DIM, DIM>;
786 pip.push_back(
787 new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
788
790}
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66

◆ 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 793 of file HenckyOps.hpp.

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

◆ 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
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 {
75 eig(i) = eig_c(i);
76 eigen_vec(i, j) = eigen_vec_c(i, j);
77 }
78};
auto is_eq(const double &a, const double &b)

Variable Documentation

◆ d_f

auto HenckyOps::d_f = [](double v) { return 0.5 / v; }

Definition at line 16 of file HenckyOps.hpp.

◆ dd_f

auto HenckyOps::dd_f = [](double v) { return -0.5 / (v * v); }

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); }

Definition at line 15 of file HenckyOps.hpp.