8#ifndef __HENKY_OPS_HPP__
9#define __HENKY_OPS_HPP__
13constexpr double eps = std::numeric_limits<float>::epsilon();
15auto f = [](
double v) {
return 0.5 * std::log(
v); };
16auto d_f = [](
double v) {
return 0.5 /
v; };
17auto dd_f = [](
double v) {
return -0.5 / (
v *
v); };
20 static inline auto check(
const double &a,
const double &b) {
28inline auto is_eq(
const double &a,
const double &b) {
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));
45 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
47 int i = 0,
j = 1,
k = 2;
49 if (
is_eq(eig(0), eig(1))) {
53 }
else if (
is_eq(eig(0), eig(2))) {
57 }
else if (
is_eq(eig(1), eig(2))) {
64 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
66 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
68 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
76 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
80struct CommonData :
public boost::enable_shared_from_this<CommonData> {
95 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
100 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
105 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matLogC);
109 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matTangent);
113template <
int DIM, IntegrationType I,
typename DomainEleOp>
116template <
int DIM, IntegrationType I,
typename DomainEleOp>
119template <
int DIM, IntegrationType I,
typename DomainEleOp>
122template <
int DIM, IntegrationType I,
typename DomainEleOp>
125template <
int DIM, IntegrationType I,
typename DomainEleOp>
128template <
int DIM, IntegrationType I,
typename DomainEleOp>
131template <
int DIM, IntegrationType I,
typename DomainEleOp>
134template <
int DIM,
typename DomainEleOp>
138 boost::shared_ptr<CommonData> common_data)
140 commonDataPtr(common_data) {
156 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
158 commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts,
false);
159 commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts,
false);
160 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
161 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
163 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
173 for (
int ii = 0; ii != DIM; ii++)
174 for (
int jj = 0; jj != DIM; jj++)
175 eigen_vec(ii, jj) = C(ii, jj);
177 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
178 for (
auto ii = 0; ii != DIM; ++ii)
179 eig(ii) = std::max(std::numeric_limits<double>::epsilon(), eig(ii));
182 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
183 if constexpr (DIM == 3) {
185 sort_eigen_vals<DIM>(eig, eigen_vec);
189 t_eig_val(
i) = eig(
i);
190 t_eig_vec(
i,
j) = eigen_vec(
i,
j);
204template <
int DIM,
typename DomainEleOp>
208 boost::shared_ptr<CommonData> common_data)
210 commonDataPtr(common_data) {
224 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
225 commonDataPtr->matLogC.resize(
size_symm, nb_gauss_pts,
false);
227 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
228 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
230 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
232 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
235 auto nb_uniq = get_uniq_nb<DIM>(&(t_eig_val(0)));
239 eig(
i) = t_eig_val(
i);
240 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
242 t_logC(
i,
j) = logC(
i,
j);
256template <
int DIM,
typename DomainEleOp>
260 boost::shared_ptr<CommonData> common_data)
262 commonDataPtr(common_data) {
276 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
278 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
279 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
280 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
282 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
286 eig(
i) = t_eig_val(
i);
287 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
290 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
292 dlogC_dC(
i,
j,
k,
l) *= 2;
294 t_logC_dC(
i,
j,
k,
l) = dlogC_dC(
i,
j,
k,
l);
308template <
int DIM,
typename DomainEleOp>
313 boost::shared_ptr<CommonData> common_data)
315 commonDataPtr(common_data) {
332 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
333 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
334 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
335 commonDataPtr->matHenckyStress.resize(
size_symm, nb_gauss_pts,
false);
336 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
338 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
339 t_T(
i,
j) = t_D(
i,
j,
k,
l) * t_logC(
k,
l);
352template <
int DIM,
typename DomainEleOp>
357 boost::shared_ptr<CommonData> common_data,
358 boost::shared_ptr<MatrixDouble> mat_D_ptr,
359 const double scale = 1)
361 scaleStress(
scale), matDPtr(mat_D_ptr) {
365 matLogCPlastic = commonDataPtr->matLogCPlastic;
380 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
381 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
382 auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
383 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
384 commonDataPtr->matHenckyStress.resize(
size_symm, nb_gauss_pts,
false);
385 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
387 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
388 t_T(
i,
j) = t_D(
i,
j,
k,
l) * (t_logC(
k,
l) - t_logCPlastic(
k,
l));
389 t_T(
i,
j) /= scaleStress;
406template <
int DIM,
typename DomainEleOp>
411 boost::shared_ptr<CommonData> common_data)
413 commonDataPtr(common_data) {
430 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
431 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
432 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
433 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
434 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts,
false);
435 commonDataPtr->matSecondPiolaStress.resize(
size_symm, nb_gauss_pts,
false);
436 auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
437 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
439 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
440 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
442 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
445 t_S(
k,
l) = t_T(
i,
j) * t_logC_dC(
i,
j,
k,
l);
446 t_P(
i,
l) = t_F(
i,
k) * t_S(
k,
l);
464template <
int DIM,
typename DomainEleOp>
467 boost::shared_ptr<CommonData> common_data,
468 boost::shared_ptr<MatrixDouble> mat_D_ptr =
nullptr)
470 commonDataPtr(common_data) {
476 matDPtr = commonDataPtr->matDPtr;
494 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
495 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
497 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
499 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*matDPtr);
500 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
501 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
502 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
504 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
505 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
506 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
508 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
516 eig(
i) = t_eig_val(
i);
517 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
521 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
531 P_D_P_plus_TL(
i,
j,
k,
l) =
533 (t_logC_dC(
i,
j,
o,
p) * t_D(
o,
p,
m,
n)) * t_logC_dC(
m,
n,
k,
l);
534 P_D_P_plus_TL(
i,
j,
k,
l) *= 0.5;
537 t_F(
i,
k) * (P_D_P_plus_TL(
k,
j,
o,
p) * dC_dF(
o,
p,
m,
n));
559 template <
int DIM, IntegrationType I>
562 template <
int DIM, IntegrationType I>
565 template <
int DIM, IntegrationType I>
568 template <
int DIM, IntegrationType I>
572 template <
int DIM, IntegrationType I>
576 template <
int DIM, IntegrationType I>
580 template <
int DIM, IntegrationType I>
587 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
588 std::string
field_name, std::string block_name,
589 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
594 OpMatBlocks(std::string
field_name, boost::shared_ptr<MatrixDouble>
m,
597 std::vector<const CubitMeshSets *> meshset_vec_ptr,
602 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]),
false);
604 "Can not get data from block");
607 MoFEMErrorCode doWork(
int side,
EntityType type,
608 EntitiesFieldData::EntData &data) {
611 for (
auto &b : blockData) {
613 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
614 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
615 b.shearModulusG * scaleYoungModulus);
620 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
621 shearModulusGDefault * scaleYoungModulus);
626 boost::shared_ptr<MatrixDouble> matDPtr;
627 const double scaleYoungModulus;
631 double shearModulusG;
635 double bulkModulusKDefault;
636 double shearModulusGDefault;
637 std::vector<BlockData> blockData;
641 std::vector<const CubitMeshSets *> meshset_vec_ptr,
645 for (
auto m : meshset_vec_ptr) {
647 std::vector<double> block_data;
648 CHKERR m->getAttributes(block_data);
649 if (block_data.size() != 2) {
651 "Expected that block has two attribute");
653 auto get_block_ents = [&]() {
656 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
675 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
679 auto set_material_stiffness = [&]() {
685 double A = (DIM == 2)
689 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
696 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
698 set_material_stiffness();
706 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"hencky_",
"",
"none");
707 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
709 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"", nu, &nu,
711 ierr = PetscOptionsEnd();
715 pip.push_back(
new OpMatBlocks(
719 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
721 (boost::format(
"%s(.*)") % block_name).str()
731template <
int DIM, IntegrationType I,
typename DomainEleOp>
734 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
735 std::string
field_name, std::string block_name, Sev sev,
double scale = 1) {
737 auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
738 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
739 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
742 common_ptr->matDPtr, sev,
scale),
749 pip.push_back(
new typename H::template OpCalculateEigenVals<DIM, I>(
752 new typename H::template OpCalculateLogC<DIM, I>(
field_name, common_ptr));
753 pip.push_back(
new typename H::template OpCalculateLogC_dC<DIM, I>(
755 pip.push_back(
new typename H::template OpCalculateHenckyStress<DIM, I>(
757 pip.push_back(
new typename H::template OpCalculatePiolaStress<DIM, I>(
763template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
766 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
767 std::string
field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
771 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
773 using OpInternalForcePiola =
776 new OpInternalForcePiola(
"U", common_ptr->getMatFirstPiolaStress()));
781template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
784 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
785 std::string
field_name, std::string block_name, Sev sev,
double scale = 1) {
788 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
796template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
799 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
800 std::string
field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
804 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
806 using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
810 new typename H::template OpHenckyTangent<DIM, I>(
field_name, common_ptr));
817template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
820 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
821 std::string
field_name, std::string block_name, Sev sev,
double scale = 1) {
824 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
static PetscErrorCode ierr
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto get_uniq_nb(double *ptr)
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
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)
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)
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)
auto is_eq(const double &a, const double &b)
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)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
double young_modulus
Young modulus.
constexpr auto field_name
auto getMatHenckyStress()
auto getMatFirstPiolaStress()
boost::shared_ptr< MatrixDouble > matLogCPlastic
MatrixDouble matFirstPiolaStress
boost::shared_ptr< MatrixDouble > matDPtr
MatrixDouble matSecondPiolaStress
MatrixDouble matHenckyStress
boost::shared_ptr< MatrixDouble > matGradPtr
OpCalculateEigenValsImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHenckyPlasticStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr, const double scale=1)
boost::shared_ptr< MatrixDouble > matLogCPlastic
boost::shared_ptr< MatrixDouble > matDPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHenckyStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateLogC_dCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateLogCImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePiolaStressImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > matDPtr
OpHenckyTangentImpl(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
boost::shared_ptr< CommonData > commonDataPtr
static auto check(const double &a, const double &b)
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.