21 boost::shared_ptr<HMHHencky> hencky_ptr)
26 "Can not get data from block");
35 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
41 b.bulkModulusK, b.shearModulusG);
73 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
74 boost::shared_ptr<PhysicalEquations> physics_ptr) {
76 data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
82 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
83 const double alpha_u);
98 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
99 const double alpha_u) {
106 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
120 std::string row_field, std::string col_field,
121 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
141 boost::shared_ptr<double> total_energy_ptr);
151 boost::shared_ptr<double> total_energy_ptr) {
157 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
158 boost::shared_ptr<HMHHencky> hencky_ptr);
162 boost::shared_ptr<DataAtIntegrationPts>
168 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
169 boost::shared_ptr<PhysicalEquations> physics_ptr) {
171 data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
176 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"hencky_",
"",
"none");
178 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
180 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"",
nu, &
nu,
183 ierr = PetscOptionsEnd();
194 (boost::format(
"%s(.*)") %
"MAT_ELASTIC").str()
206 for (
auto m : meshset_vec_ptr) {
208 std::vector<double> block_data;
209 CHKERR m->getAttributes(block_data);
210 if (block_data.size() != 2) {
212 "Expected that block has two attribute");
214 auto get_block_ents = [&]() {
235 boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
236 boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
240 auto set_material_stiffness = [&]() {
251 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
252 &*(mat_D_ptr->data().begin()));
253 auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
254 &*mat_axiator_D_ptr->data().begin());
255 auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
256 &*mat_deviator_D_ptr->data().begin());
257 t_axiator_D(
i,
j,
k,
l) =
A *
260 t_deviator_D(
i,
j,
k,
l) =
262 t_D(
i,
j,
k,
l) = t_axiator_D(
i,
j,
k,
l) + t_deviator_D(
i,
j,
k,
l);
269 set_material_stiffness();
289 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u)
294 "get polyconvex option failed");
300 CHKERR integratePolyconvexHencky(data);
302 CHKERR integrateHencky(data);
314 int nb_integration_pts = data.
getN().size1();
315 auto v = getVolume();
316 auto t_w = getFTensor0IntegrationWeight();
317 auto t_approx_P_adjont_log_du =
318 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
319 auto t_log_stretch_h1 =
320 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
322 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
324 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
330 auto get_ftensor2 = [](
auto &
v) {
332 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
335 int nb_base_functions = data.
getN().size2();
337 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
339 auto t_nf = get_ftensor2(nF);
343 t_D(
i,
j,
k,
l) * (t_log_stretch_h1(
k,
l) + alphaU * t_dot_log_u(
k,
l));
346 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
349 for (; bb != nb_dofs / 6; ++bb) {
350 t_nf(
L) += t_row_base_fun * t_residual(
L);
354 for (; bb != nb_base_functions; ++bb)
358 ++t_approx_P_adjont_log_du;
373 int nb_integration_pts = data.
getN().size1();
374 auto v = getVolume();
375 auto t_w = getFTensor0IntegrationWeight();
376 auto t_approx_P_adjont_log_du =
377 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
378 auto t_log_stretch_h1 =
379 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
381 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
383 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
389 auto get_ftensor2 = [](
auto &
v) {
391 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
394 constexpr
double nohat_k = 1. / 4;
395 constexpr
double hat_k = 1. / 8;
396 double mu = dataAtPts->mu;
397 double lambda = dataAtPts->lambda;
399 constexpr
double third = boost::math::constants::third<double>();
403 int nb_base_functions = data.
getN().size2();
405 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
407 auto t_nf = get_ftensor2(nF);
409 double log_det = t_log_stretch_h1(
i,
i);
410 double log_det2 = log_det * log_det;
413 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
416 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
417 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
420 A * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j))
428 alphaU * t_D(
i,
j,
k,
l) * t_dot_log_u(
k,
l);
432 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
435 for (; bb != nb_dofs /
size_symm; ++bb) {
436 t_nf(
L) += t_row_base_fun * t_residual(
L);
440 for (; bb != nb_base_functions; ++bb)
444 ++t_approx_P_adjont_log_du;
452 std::string row_field, std::string col_field,
453 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha)
460 "get polyconvex option failed");
469 CHKERR integratePolyconvexHencky(row_data, col_data);
471 CHKERR integrateHencky(row_data, col_data);
486 int nb_integration_pts = row_data.
getN().size1();
487 int row_nb_dofs = row_data.
getIndices().size();
488 int col_nb_dofs = col_data.
getIndices().size();
494 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
495 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
497 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
498 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
500 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
501 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
503 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
504 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
506 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
507 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
509 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
510 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
521 auto v = getVolume();
522 auto t_w = getFTensor0IntegrationWeight();
524 auto t_approx_P_adjont_dstretch =
525 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
526 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
527 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
528 auto &nbUniq = dataAtPts->nbUniq;
530 int row_nb_base_functions = row_data.
getN().size2();
533 auto get_dP = [&]() {
535 auto ts_a = getTSa();
537 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
539 t_dP_tmp(
L,
J) = -(1 + alphaU * ts_a) *
541 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
544 auto t_approx_P_adjont_dstretch =
545 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
546 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
547 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
548 auto &nbUniq = dataAtPts->nbUniq;
550 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
551 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
556 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
557 t_approx_P_adjont_dstretch(
j,
i));
562 t_dP(
L,
J) = t_L(
i,
j,
L) *
563 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
569 ++t_approx_P_adjont_dstretch;
574 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
575 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
576 t_dP(
L,
J) = t_dP_tmp(
L,
J);
581 return getFTensor2FromMat<size_symm, size_symm>(dP);
584 auto t_dP = get_dP();
585 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
589 for (; rr != row_nb_dofs / 6; ++rr) {
591 auto t_m = get_ftensor2(
K, 6 * rr, 0);
592 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
593 const double b =
a * t_row_base_fun * t_col_base_fun;
594 t_m(
L,
J) += b * t_dP(
L,
J);
601 for (; rr != row_nb_base_functions; ++rr) {
620 int nb_integration_pts = row_data.
getN().size1();
621 int row_nb_dofs = row_data.
getIndices().size();
622 int col_nb_dofs = col_data.
getIndices().size();
628 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
629 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
631 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
632 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
634 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
635 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
637 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
638 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
640 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
641 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
643 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
644 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
655 auto v = getVolume();
656 auto t_w = getFTensor0IntegrationWeight();
658 int row_nb_base_functions = row_data.
getN().size2();
661 auto get_dP = [&]() {
663 auto ts_a = getTSa();
665 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
667 constexpr
double nohat_k = 1. / 4;
668 constexpr
double hat_k = 1. / 8;
669 double mu = dataAtPts->mu;
670 double lambda = dataAtPts->lambda;
672 constexpr
double third = boost::math::constants::third<double>();
676 auto t_approx_P_adjont_dstretch =
677 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
678 auto t_log_stretch_h1 =
679 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
680 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
681 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
682 auto &nbUniq = dataAtPts->nbUniq;
684 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
685 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
687 double log_det = t_log_stretch_h1(
i,
i);
688 double log_det2 = log_det * log_det;
691 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
693 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
694 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
698 (
A * 2 * nohat_k) * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j));
699 t_B_diff(
i,
j) = (B * 2 * hat_k) * log_det *
t_kd(
i,
j) +
703 t_A_diff(
i,
j) * (t_dev(
m,
n) * t_diff_deviator(
m,
n,
k,
l))
707 A * t_diff_deviator(
m,
n,
i,
j) * t_diff_deviator(
m,
n,
k,
l)
713 t_dP(
L,
J) = -t_L(
i,
j,
L) *
720 (alphaU * ts_a) * (t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)
729 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
730 t_approx_P_adjont_dstretch(
j,
i));
735 t_dP(
L,
J) += t_L(
i,
j,
L) *
736 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
742 ++t_approx_P_adjont_dstretch;
748 return getFTensor2FromMat<size_symm, size_symm>(dP);
751 auto t_dP = get_dP();
752 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
756 for (; rr != row_nb_dofs / 6; ++rr) {
758 auto t_m = get_ftensor2(
K, 6 * rr, 0);
759 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
760 const double b =
a * t_row_base_fun * t_col_base_fun;
761 t_m(
L,
J) += b * t_dP(
L,
J);
768 for (; rr != row_nb_base_functions; ++rr) {
779 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
780 boost::shared_ptr<double> total_energy_ptr)
782 totalEnergyPtr(total_energy_ptr) {}
793 int nb_integration_pts = getGaussPts().size2();
795 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
796 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
798 dataAtPts->energyAtPts.resize(nb_integration_pts,
false);
801 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
803 t_energy = 0.5 * (t_log_u(
i,
j) * (t_D(
i,
j,
k,
l) * t_log_u(
k,
l)));
809 if (totalEnergyPtr) {
810 auto t_w = getFTensor0IntegrationWeight();
812 double loc_energy = 0;
813 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
814 loc_energy += t_energy * t_w;
818 *totalEnergyPtr += getMeasure() * loc_energy;
825 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
826 boost::shared_ptr<HMHHencky> hencky_ptr)
828 henckyPtr(hencky_ptr) {}
842 auto nb_integration_pts = dataAtPts->approxPAtPts.size2();
844 if (nb_integration_pts != getGaussPts().size2()) {
846 "inconsistent number of integration points");
852 for (
auto &b : henckyPtr->blockData) {
854 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
855 dataAtPts->mu = b.bulkModulusK;
856 dataAtPts->lambda = b.bulkModulusK;
857 CHKERR henckyPtr->getMatDPtr(
858 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
859 dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
864 const auto E = henckyPtr->E;
865 const auto nu = henckyPtr->nu;
873 CHKERR henckyPtr->getMatDPtr(
874 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
879 auto get_invert_D = [&]() {
882 noalias(dataAtPts->matInvD) = dataAtPts->matD;
890 dataAtPts->logStretchTensorAtPts.resize(
size_symm, nb_integration_pts,
false);
893 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
895 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
896 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
897 &*dataAtPts->matInvD.data().begin());
901 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
902 t_log_u(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_approx_P(
k,
l);