21 boost::shared_ptr<HMHHencky> hencky_ptr)
26 "Can not get data from block");
34 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
40 b.bulkModulusK, b.shearModulusG);
71 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
72 boost::shared_ptr<PhysicalEquations> physics_ptr) {
74 data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
80 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
81 const double alpha_u);
96 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
97 const double alpha_u) {
104 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
118 std::string row_field, std::string col_field,
119 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
139 boost::shared_ptr<double> total_energy_ptr);
149 boost::shared_ptr<double> total_energy_ptr) {
155 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
156 boost::shared_ptr<MatrixDouble> strain_ptr,
157 boost::shared_ptr<MatrixDouble> stress_ptr,
158 boost::shared_ptr<HMHHencky> hencky_ptr);
162 boost::shared_ptr<DataAtIntegrationPts>
170 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
171 boost::shared_ptr<PhysicalEquations> physics_ptr) {
173 data_ptr, data_ptr->getLogStretchTensorAtPts(),
174 data_ptr->getApproxPAtPts(),
175 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
179 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
180 boost::shared_ptr<PhysicalEquations> physics_ptr) {
182 data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
183 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
188 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"hencky_",
"",
"none");
190 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
192 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"",
nu, &
nu,
196 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Hencky: E = " <<
E <<
" nu = " <<
nu;
198 ierr = PetscOptionsEnd();
209 (boost::format(
"%s(.*)") %
"MAT_ELASTIC").str()
221 for (
auto m : meshset_vec_ptr) {
223 std::vector<double> block_data;
224 CHKERR m->getAttributes(block_data);
225 if (block_data.size() != 2) {
227 "Expected that block has two attribute");
229 auto get_block_ents = [&]() {
250 boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
251 boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
255 auto set_material_stiffness = [&]() {
261 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
262 &*(mat_D_ptr->data().begin()));
263 auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
264 &*mat_axiator_D_ptr->data().begin());
265 auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
266 &*mat_deviator_D_ptr->data().begin());
269 t_deviator_D(
i,
j,
k,
l) =
271 t_D(
i,
j,
k,
l) = t_axiator_D(
i,
j,
k,
l) + t_deviator_D(
i,
j,
k,
l);
278 set_material_stiffness();
287 auto set_material_compilance = [&]() {
296 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
297 &*(mat_inv_D_ptr->data().begin()));
298 t_inv_D(
i,
j,
k,
l) =
304 set_material_compilance();
324 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u)
329 "get polyconvex option failed");
335 CHKERR integratePolyconvexHencky(data);
337 CHKERR integrateHencky(data);
349 int nb_integration_pts = data.
getN().size1();
350 auto v = getVolume();
351 auto t_w = getFTensor0IntegrationWeight();
352 auto t_approx_P_adjoint_log_du =
353 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
354 auto t_log_stretch_h1 =
355 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
357 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
359 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
365 auto get_ftensor2 = [](
auto &
v) {
367 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
370 int nb_base_functions = data.
getN().size2();
372 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
374 auto t_nf = get_ftensor2(nF);
378 t_D(
i,
j,
k,
l) * (t_log_stretch_h1(
k,
l) + alphaU * t_dot_log_u(
k,
l));
381 a * (t_approx_P_adjoint_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
384 for (; bb != nb_dofs / 6; ++bb) {
385 t_nf(
L) -= t_row_base_fun * t_residual(
L);
389 for (; bb != nb_base_functions; ++bb)
393 ++t_approx_P_adjoint_log_du;
408 int nb_integration_pts = data.
getN().size1();
409 auto v = getVolume();
410 auto t_w = getFTensor0IntegrationWeight();
411 auto t_approx_P_adjoint_log_du =
412 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
413 auto t_log_stretch_h1 =
414 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
416 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
418 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
424 auto get_ftensor2 = [](
auto &
v) {
426 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
429 constexpr
double nohat_k = 1. / 4;
430 constexpr
double hat_k = 1. / 8;
431 double mu = dataAtPts->mu;
432 double lambda = dataAtPts->lambda;
434 constexpr
double third = boost::math::constants::third<double>();
438 int nb_base_functions = data.
getN().size2();
440 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
442 auto t_nf = get_ftensor2(nF);
444 double log_det = t_log_stretch_h1(
i,
i);
445 double log_det2 = log_det * log_det;
448 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
451 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
452 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
455 A * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j))
463 alphaU * t_D(
i,
j,
k,
l) * t_dot_log_u(
k,
l);
467 a * (t_approx_P_adjoint_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
470 for (; bb != nb_dofs /
size_symm; ++bb) {
471 t_nf(
L) -= t_row_base_fun * t_residual(
L);
475 for (; bb != nb_base_functions; ++bb)
479 ++t_approx_P_adjoint_log_du;
487 std::string row_field, std::string col_field,
488 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha)
495 "get polyconvex option failed");
504 CHKERR integratePolyconvexHencky(row_data, col_data);
506 CHKERR integrateHencky(row_data, col_data);
521 int nb_integration_pts = row_data.
getN().size1();
522 int row_nb_dofs = row_data.
getIndices().size();
523 int col_nb_dofs = col_data.
getIndices().size();
529 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
530 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
532 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
533 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
535 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
536 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
538 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
539 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
541 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
542 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
544 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
545 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
556 auto v = getVolume();
557 auto t_w = getFTensor0IntegrationWeight();
559 auto t_approx_P_adjoint_dstretch =
560 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
561 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
562 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
563 auto &nbUniq = dataAtPts->nbUniq;
565 int row_nb_base_functions = row_data.
getN().size2();
568 auto get_dP = [&]() {
570 auto ts_a = getTSa();
572 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
574 t_dP_tmp(
L,
J) = -(1 + alphaU * ts_a) *
576 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
580 auto t_approx_P_adjoint_dstretch =
581 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
582 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
583 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
584 auto &nbUniq = dataAtPts->nbUniq;
586 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
587 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
592 t_sym(
i,
j) = (t_approx_P_adjoint_dstretch(
i,
j) ||
593 t_approx_P_adjoint_dstretch(
j,
i));
598 t_dP(
L,
J) = t_L(
i,
j,
L) *
599 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
605 ++t_approx_P_adjoint_dstretch;
610 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
611 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
612 t_dP(
L,
J) = t_dP_tmp(
L,
J);
617 return getFTensor2FromMat<size_symm, size_symm>(dP);
620 auto t_dP = get_dP();
621 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
625 for (; rr != row_nb_dofs / 6; ++rr) {
627 auto t_m = get_ftensor2(
K, 6 * rr, 0);
628 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
629 const double b =
a * t_row_base_fun * t_col_base_fun;
630 t_m(
L,
J) -= b * t_dP(
L,
J);
637 for (; rr != row_nb_base_functions; ++rr) {
656 int nb_integration_pts = row_data.
getN().size1();
657 int row_nb_dofs = row_data.
getIndices().size();
658 int col_nb_dofs = col_data.
getIndices().size();
664 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
665 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
667 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
668 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
670 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
671 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
673 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
674 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
676 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
677 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
679 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
680 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
691 auto v = getVolume();
692 auto t_w = getFTensor0IntegrationWeight();
694 int row_nb_base_functions = row_data.
getN().size2();
697 auto get_dP = [&]() {
699 auto ts_a = getTSa();
701 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
703 constexpr
double nohat_k = 1. / 4;
704 constexpr
double hat_k = 1. / 8;
705 double mu = dataAtPts->mu;
706 double lambda = dataAtPts->lambda;
708 constexpr
double third = boost::math::constants::third<double>();
712 auto t_approx_P_adjoint_dstretch =
713 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
714 auto t_log_stretch_h1 =
715 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
716 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
717 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
718 auto &nbUniq = dataAtPts->nbUniq;
720 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
721 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
723 double log_det = t_log_stretch_h1(
i,
i);
724 double log_det2 = log_det * log_det;
727 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
729 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
730 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
734 (
A * 2 * nohat_k) * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j));
735 t_B_diff(
i,
j) = (B * 2 * hat_k) * log_det *
t_kd(
i,
j) +
739 t_A_diff(
i,
j) * (t_dev(
m,
n) * t_diff_deviator(
m,
n,
k,
l))
743 A * t_diff_deviator(
m,
n,
i,
j) * t_diff_deviator(
m,
n,
k,
l)
749 t_dP(
L,
J) = -t_L(
i,
j,
L) *
756 (alphaU * ts_a) * (t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)
766 t_sym(
i,
j) = (t_approx_P_adjoint_dstretch(
i,
j) ||
767 t_approx_P_adjoint_dstretch(
j,
i));
772 t_dP(
L,
J) += t_L(
i,
j,
L) *
773 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
779 ++t_approx_P_adjoint_dstretch;
785 return getFTensor2FromMat<size_symm, size_symm>(dP);
788 auto t_dP = get_dP();
789 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
793 for (; rr != row_nb_dofs / 6; ++rr) {
795 auto t_m = get_ftensor2(
K, 6 * rr, 0);
796 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
797 const double b =
a * t_row_base_fun * t_col_base_fun;
798 t_m(
L,
J) -= b * t_dP(
L,
J);
805 for (; rr != row_nb_base_functions; ++rr) {
816 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
817 boost::shared_ptr<double> total_energy_ptr)
819 totalEnergyPtr(total_energy_ptr) {}
830 int nb_integration_pts = getGaussPts().size2();
832 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
833 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
835 dataAtPts->energyAtPts.resize(nb_integration_pts,
false);
838 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
840 t_energy = 0.5 * (t_log_u(
i,
j) * (t_D(
i,
j,
k,
l) * t_log_u(
k,
l)));
846 if (totalEnergyPtr) {
847 auto t_w = getFTensor0IntegrationWeight();
849 double loc_energy = 0;
850 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
851 loc_energy += t_energy * t_w;
855 *totalEnergyPtr += getMeasure() * loc_energy;
862 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
863 boost::shared_ptr<MatrixDouble> strain_ptr,
864 boost::shared_ptr<MatrixDouble> stress_ptr,
865 boost::shared_ptr<HMHHencky> hencky_ptr)
867 strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {}
881 auto nb_integration_pts = stressPtr->size2();
883 if (nb_integration_pts != getGaussPts().size2()) {
885 "inconsistent number of integration points");
891 for (
auto &b : henckyPtr->blockData) {
893 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
894 dataAtPts->mu = b.shearModulusG;
895 dataAtPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
897 CHKERR henckyPtr->getMatDPtr(
898 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
899 dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
901 CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(),
902 b.bulkModulusK, b.shearModulusG);
907 const auto E = henckyPtr->E;
908 const auto nu = henckyPtr->nu;
916 CHKERR henckyPtr->getMatDPtr(
917 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
926 strainPtr->resize(
size_symm, nb_integration_pts,
false);
927 auto t_strain = getFTensor2SymmetricFromMat<3>(*strainPtr);
928 auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressPtr);
929 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
930 &*dataAtPts->matInvD.data().begin());
932 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
933 &*dataAtPts->matD.data().begin());
936 const double lambda = dataAtPts->lambda;
937 const double mu = dataAtPts->mu;
942 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
943 t_strain(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_stress(
k,
l);
947 t_stress_symm_debug(
i,
j) = (t_stress(
i,
j) || t_stress(
j,
i)) / 2;
949 t_stress_symm_debug_diff(
i,
j) =
950 t_D(
i,
j,
k,
l) * t_strain(
k,
l) - t_stress_symm_debug(
i,
j);
952 t_stress_symm_debug_diff(
i,
j) * t_stress_symm_debug_diff(
i,
j);
953 double nrm0 = t_stress_symm_debug(
i,
j) * t_stress_symm_debug(
i,
j) +
954 std::numeric_limits<double>::epsilon();
955 constexpr
double eps = 1e-10;
956 if (std::fabs(std::sqrt(nrm / nrm0)) >
eps) {
958 <<
"Stress symmetry check failed: " << std::endl
959 << t_stress_symm_debug_diff << std::endl
962 "Norm is too big: " + std::to_string(nrm / nrm0));