8 #ifndef __HENCKY_OPS_HPP__
9 #define __HENCKY_OPS_HPP__
13 constexpr
double eps = std::numeric_limits<float>::epsilon();
15 auto f = [](
double v) {
return 0.5 * std::log(
v); };
16 auto d_f = [](
double v) {
return 0.5 /
v; };
17 auto dd_f = [](
double v) {
return -0.5 / (
v *
v); };
20 static inline auto check(
const double &
a,
const double &b) {
28 inline 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);
80 struct 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);
113 template <
int DIM, IntegrationType I,
typename DomainEleOp>
116 template <
int DIM, IntegrationType I,
typename DomainEleOp>
119 template <
int DIM, IntegrationType I,
typename DomainEleOp>
122 template <
int DIM, IntegrationType I,
typename DomainEleOp,
int S>
125 template <
int DIM, IntegrationType I,
typename DomainEleOp,
int S>
128 template <
int DIM, IntegrationType I,
typename DomainEleOp,
int S>
131 template <
int DIM, IntegrationType I,
typename DomainEleOp,
int S>
134 template <
int DIM,
typename DomainEleOp>
138 boost::shared_ptr<CommonData> common_data)
140 commonDataPtr(common_data) {
141 std::fill(&DomainEleOp::doEntities[MBEDGE],
142 &DomainEleOp::doEntities[MBMAXTYPE],
false);
154 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
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);
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);
204 template <
int DIM,
typename DomainEleOp>
208 boost::shared_ptr<CommonData> common_data)
210 commonDataPtr(common_data) {
211 std::fill(&DomainEleOp::doEntities[MBEDGE],
212 &DomainEleOp::doEntities[MBMAXTYPE],
false);
222 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
223 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
224 commonDataPtr->matLogC.resize(
size_symm, nb_gauss_pts,
false);
226 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
227 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
229 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
231 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
245 template <
int DIM,
typename DomainEleOp>
249 boost::shared_ptr<CommonData> common_data)
251 commonDataPtr(common_data) {
252 std::fill(&DomainEleOp::doEntities[MBEDGE],
253 &DomainEleOp::doEntities[MBMAXTYPE],
false);
264 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
265 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
267 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
268 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
269 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
271 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
273 auto nb_uniq = get_uniq_nb<DIM>(&t_eig_val(0));
274 t_logC_dC(
i,
j,
k,
l) =
276 nb_uniq)(
i,
j,
k,
l);
290 template <
int DIM,
typename DomainEleOp,
int S>
295 boost::shared_ptr<CommonData> common_data)
297 commonDataPtr(common_data) {
298 std::fill(&DomainEleOp::doEntities[MBEDGE],
299 &DomainEleOp::doEntities[MBMAXTYPE],
false);
311 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
312 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
313 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
314 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
315 commonDataPtr->matHenckyStress.resize(
size_symm, nb_gauss_pts,
false);
316 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
318 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
319 t_T(
i,
j) = t_D(
i,
j,
k,
l) * t_logC(
k,
l);
332 template <
int DIM,
typename DomainEleOp,
int S>
337 boost::shared_ptr<CommonData> common_data,
338 boost::shared_ptr<MatrixDouble> mat_D_ptr,
339 const double scale = 1)
341 scaleStress(
scale), matDPtr(mat_D_ptr) {
342 std::fill(&DomainEleOp::doEntities[MBEDGE],
343 &DomainEleOp::doEntities[MBMAXTYPE],
false);
345 matLogCPlastic = commonDataPtr->matLogCPlastic;
357 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
358 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
359 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
360 auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
361 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
362 commonDataPtr->matHenckyStress.resize(
size_symm, nb_gauss_pts,
false);
363 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
365 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
366 t_T(
i,
j) = t_D(
i,
j,
k,
l) * (t_logC(
k,
l) - t_logCPlastic(
k,
l));
367 t_T(
i,
j) /= scaleStress;
384 template <
int DIM,
typename DomainEleOp,
int S>
389 boost::shared_ptr<CommonData> common_data)
391 commonDataPtr(common_data) {
392 std::fill(&DomainEleOp::doEntities[MBEDGE],
393 &DomainEleOp::doEntities[MBMAXTYPE],
false);
407 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
408 #ifdef HENCKY_SMALL_STRAIN
409 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
411 auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
412 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
413 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
414 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts,
false);
415 commonDataPtr->matSecondPiolaStress.resize(
size_symm, nb_gauss_pts,
false);
416 auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
417 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
419 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
420 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
422 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
424 #ifdef HENCKY_SMALL_STRAIN
425 t_P(
i,
j) = t_D(
i,
j,
k,
l) * t_grad(
k,
l);
429 t_S(
k,
l) = t_T(
i,
j) * t_logC_dC(
i,
j,
k,
l);
430 t_P(
i,
l) = t_F(
i,
k) * t_S(
k,
l);
439 #ifdef HENCKY_SMALL_STRAIN
451 template <
int DIM,
typename DomainEleOp,
int S>
454 boost::shared_ptr<CommonData> common_data,
455 boost::shared_ptr<MatrixDouble> mat_D_ptr =
nullptr)
457 commonDataPtr(common_data) {
458 std::fill(&DomainEleOp::doEntities[MBEDGE],
459 &DomainEleOp::doEntities[MBMAXTYPE],
false);
463 matDPtr = commonDataPtr->matDPtr;
480 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
481 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
483 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
485 auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
486 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
487 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
488 auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
490 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
491 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
492 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
494 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
496 #ifdef HENCKY_SMALL_STRAIN
497 dP_dF(
i,
j,
k,
l) = t_D(
i,
j,
k,
l);
504 auto nb_uniq = get_uniq_nb<DIM>(&t_eig_val(0));
512 P_D_P_plus_TL(
i,
j,
k,
l) =
514 (t_logC_dC(
i,
j, o, p) * t_D(o, p,
m,
n)) * t_logC_dC(
m,
n,
k,
l);
515 P_D_P_plus_TL(
i,
j,
k,
l) *= 0.5;
518 t_F(
i,
k) * (P_D_P_plus_TL(
k,
j, o, p) * dC_dF(o, p,
m,
n));
542 template <
int DIM, IntegrationType I>
545 template <
int DIM, IntegrationType I>
548 template <
int DIM, IntegrationType I>
551 template <
int DIM, IntegrationType I,
int S>
555 template <
int DIM, IntegrationType I,
int S = 0>
559 template <
int DIM, IntegrationType I,
int S>
563 template <
int DIM, IntegrationType I,
int S>
570 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
571 std::string block_name,
572 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev,
579 std::vector<const CubitMeshSets *> meshset_vec_ptr,
585 "Can not get data from block");
592 for (
auto &b : blockData) {
594 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
595 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
596 b.shearModulusG * scaleYoungModulus);
601 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
602 shearModulusGDefault * scaleYoungModulus);
607 boost::shared_ptr<MatrixDouble> matDPtr;
608 const double scaleYoungModulus;
612 double shearModulusG;
616 double bulkModulusKDefault;
617 double shearModulusGDefault;
618 std::vector<BlockData> blockData;
622 std::vector<const CubitMeshSets *> meshset_vec_ptr,
626 for (
auto m : meshset_vec_ptr) {
628 std::vector<double> block_data;
629 CHKERR m->getAttributes(block_data);
630 if (block_data.size() != 2) {
632 "Expected that block has two attribute");
634 auto get_block_ents = [&]() {
637 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
656 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
660 auto set_material_stiffness = [&]() {
666 double A = (DIM == 2)
670 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
677 constexpr
auto size_symm = (DIM * (DIM + 1)) / 2;
679 set_material_stiffness();
687 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"",
"none");
688 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
690 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"", nu, &nu,
692 ierr = PetscOptionsEnd();
696 pip.push_back(
new OpMatBlocks(
700 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
702 (boost::format(
"%s(.*)") % block_name).str()
712 template <
int DIM, IntegrationType I,
typename DomainEleOp>
715 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
718 auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
719 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
720 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
723 common_ptr->matDPtr, sev,
scale),
730 pip.push_back(
new typename H::template OpCalculateEigenVals<DIM, I>(
733 new typename H::template OpCalculateLogC<DIM, I>(
field_name, common_ptr));
734 pip.push_back(
new typename H::template OpCalculateLogC_dC<DIM, I>(
737 pip.push_back(
new typename H::template OpCalculateHenckyStress<DIM, I, 0>(
739 pip.push_back(
new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
745 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
748 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
749 std::string
field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
753 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
754 A>::template LinearForm<I>;
763 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
766 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
770 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
778 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
781 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
782 std::string
field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
786 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
788 using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
793 new typename H::template OpHenckyTangent<DIM, I, 0>(
field_name, common_ptr));
800 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
803 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
807 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
816 #endif // __HENCKY_OPS_HPP__