21 boost::shared_ptr<HMHHencky> hencky_ptr)
26 "Can not get data from block");
29 MoFEMErrorCode
doWork(
int side, EntityType type,
30 EntitiesFieldData::EntData &data) {
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,
105 boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
106 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
117 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
118 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
119 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
122 field_name, data_ptr, external_strain_vec_ptr, smv);
128 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
142 std::string row_field, std::string col_field,
143 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
163 boost::shared_ptr<double> total_energy_ptr);
164 MoFEMErrorCode
doWork(
int side, EntityType type,
EntData &data);
173 boost::shared_ptr<double> total_energy_ptr) {
179 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
180 boost::shared_ptr<MatrixDouble> strain_ptr,
181 boost::shared_ptr<MatrixDouble> stress_ptr,
182 boost::shared_ptr<HMHHencky> hencky_ptr);
183 MoFEMErrorCode
doWork(
int side, EntityType type,
EntData &data);
186 boost::shared_ptr<DataAtIntegrationPts>
194 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
195 boost::shared_ptr<PhysicalEquations> physics_ptr) {
197 data_ptr, data_ptr->getLogStretchTensorAtPts(),
198 data_ptr->getApproxPAtPts(),
199 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
203 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
204 boost::shared_ptr<PhysicalEquations> physics_ptr) {
206 data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
207 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
210 MoFEMErrorCode
getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
212 PetscOptionsBegin(PETSC_COMM_WORLD,
"hencky_",
"",
"none");
214 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
216 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"",
nu, &
nu,
222 <<
"Hencky: E = " <<
E <<
" nu = " <<
nu;
235 (boost::format(
"%s(.*)") %
"MAT_ELASTIC").str()
247 for (
auto m : meshset_vec_ptr) {
249 std::vector<double> block_data;
250 CHKERR m->getAttributes(block_data);
251 if (block_data.size() < 2) {
253 "Expected that block has atleast two attributes");
255 auto get_block_ents = [&]() {
275 MoFEMErrorCode
getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
276 boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
277 boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
281 auto set_material_stiffness = [&]() {
287 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
288 &*(mat_D_ptr->data().begin()));
289 auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
290 &*mat_axiator_D_ptr->data().begin());
291 auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
292 &*mat_deviator_D_ptr->data().begin());
295 t_deviator_D(
i,
j,
k,
l) =
297 t_D(
i,
j,
k,
l) = t_axiator_D(
i,
j,
k,
l) + t_deviator_D(
i,
j,
k,
l);
304 set_material_stiffness();
313 auto set_material_compilance = [&]() {
322 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
323 &*(mat_inv_D_ptr->data().begin()));
324 t_inv_D(
i,
j,
k,
l) =
330 set_material_compilance();
354 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u)
357 CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR,
"",
"-poly_convex",
359 "get polyconvex option failed");
365 CHKERR integratePolyconvexHencky(data);
367 CHKERR integrateHencky(data);
379 int nb_integration_pts = data.
getN().size1();
380 auto v = getVolume();
381 auto t_w = getFTensor0IntegrationWeight();
382 auto t_approx_P_adjoint_log_du =
383 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
384 auto t_log_stretch_h1 =
385 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
387 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
389 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
395 auto get_ftensor2 = [](
auto &
v) {
397 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
400 int nb_base_functions = data.
getN().size2();
402 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
404 auto t_nf = get_ftensor2(nF);
408 t_D(
i,
j,
k,
l) * (t_log_stretch_h1(
k,
l) + alphaU * t_dot_log_u(
k,
l));
411 a * (t_approx_P_adjoint_log_du(L) - t_L(
i,
j, L) * t_T(
i,
j));
414 for (; bb != nb_dofs / 6; ++bb) {
415 t_nf(L) -= t_row_base_fun * t_residual(L);
419 for (; bb != nb_base_functions; ++bb)
423 ++t_approx_P_adjoint_log_du;
438 int nb_integration_pts = data.
getN().size1();
439 auto v = getVolume();
440 auto t_w = getFTensor0IntegrationWeight();
441 auto t_approx_P_adjoint_log_du =
442 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
443 auto t_log_stretch_h1 =
444 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
446 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
448 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
454 auto get_ftensor2 = [](
auto &
v) {
456 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
459 constexpr double nohat_k = 1. / 4;
460 constexpr double hat_k = 1. / 8;
461 double mu = dataAtPts->mu;
462 double lambda = dataAtPts->lambda;
464 constexpr double third = boost::math::constants::third<double>();
468 int nb_base_functions = data.
getN().size2();
470 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
472 auto t_nf = get_ftensor2(nF);
474 double log_det = t_log_stretch_h1(
i,
i);
475 double log_det2 = log_det * log_det;
478 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
481 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
482 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
485 A * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j))
493 alphaU * t_D(
i,
j,
k,
l) * t_dot_log_u(
k,
l);
497 a * (t_approx_P_adjoint_log_du(L) - t_L(
i,
j, L) * t_T(
i,
j));
500 for (; bb != nb_dofs /
size_symm; ++bb) {
501 t_nf(L) -= t_row_base_fun * t_residual(L);
505 for (; bb != nb_base_functions; ++bb)
509 ++t_approx_P_adjoint_log_du;
517 std::string row_field, std::string col_field,
518 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha)
523 CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR,
"",
"-poly_convex",
525 "get polyconvex option failed");
530 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
531 boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
532 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
534 externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap(smv) {}
550 for (
auto &ext_strain_block : (*externalStrainVecPtr)) {
552 if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
555 if (scalingMethodsMap.find(ext_strain_block.blockName) !=
556 scalingMethodsMap.end()) {
558 scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
561 <<
"No scaling method found for " << ext_strain_block.blockName;
565 double external_strain_val =
scale * ext_strain_block.val;
570 int nb_integration_pts = data.
getN().size1();
571 auto v = getVolume();
572 auto t_w = getFTensor0IntegrationWeight();
580 auto get_ftensor2 = [](
auto &
v) {
582 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
585 int nb_base_functions = data.
getN().size2();
587 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
588 auto tr = 3.0 * external_strain_val;
590 auto t_nf = get_ftensor2(nF);
597 t_residual(L) =
a * (t_L(
i,
j, L) * t_T(
i,
j));
600 for (; bb != nb_dofs / 6; ++bb) {
601 t_nf(L) += t_row_base_fun * t_residual(L);
605 for (; bb != nb_base_functions; ++bb)
622 CHKERR integratePolyconvexHencky(row_data, col_data);
624 CHKERR integrateHencky(row_data, col_data);
639 int nb_integration_pts = row_data.
getN().size1();
640 int row_nb_dofs = row_data.
getIndices().size();
641 int col_nb_dofs = col_data.
getIndices().size();
643 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
647 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
648 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
650 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
651 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
653 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
654 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
656 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
657 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
659 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
660 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
662 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
663 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
674 auto v = getVolume();
675 auto t_w = getFTensor0IntegrationWeight();
677 auto t_approx_P_adjoint__dstretch =
678 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
679 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
680 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
682 int row_nb_base_functions = row_data.
getN().size2();
685 auto get_dP = [&]() {
687 auto ts_a = getTSa();
689 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
691 t_dP_tmp(L,
J) = -(1 + alphaU * ts_a) *
693 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
697 auto t_approx_P_adjoint__dstretch =
698 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
699 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
700 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
701 auto &nbUniq = dataAtPts->nbUniq;
703 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
704 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
709 t_sym(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
710 t_approx_P_adjoint__dstretch(
j,
i));
715 t_dP(L,
J) = t_L(
i,
j, L) *
716 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
722 ++t_approx_P_adjoint__dstretch;
727 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
728 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
729 t_dP(L,
J) = t_dP_tmp(L,
J);
734 return getFTensor2FromMat<size_symm, size_symm>(dP);
737 auto t_dP = get_dP();
738 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
742 for (; rr != row_nb_dofs / 6; ++rr) {
744 auto t_m = get_ftensor2(K, 6 * rr, 0);
745 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
746 const double b =
a * t_row_base_fun * t_col_base_fun;
747 t_m(L,
J) -= b * t_dP(L,
J);
754 for (; rr != row_nb_base_functions; ++rr) {
773 int nb_integration_pts = row_data.
getN().size1();
774 int row_nb_dofs = row_data.
getIndices().size();
775 int col_nb_dofs = col_data.
getIndices().size();
777 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
781 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
782 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
784 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
785 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
787 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
788 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
790 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
791 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
793 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
794 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
796 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
797 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
808 auto v = getVolume();
809 auto t_w = getFTensor0IntegrationWeight();
811 int row_nb_base_functions = row_data.
getN().size2();
814 auto get_dP = [&]() {
816 auto ts_a = getTSa();
818 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
820 constexpr double nohat_k = 1. / 4;
821 constexpr double hat_k = 1. / 8;
822 double mu = dataAtPts->mu;
823 double lambda = dataAtPts->lambda;
825 constexpr double third = boost::math::constants::third<double>();
829 auto t_approx_P_adjoint__dstretch =
830 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
831 auto t_log_stretch_h1 =
832 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
833 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
834 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
835 auto &nbUniq = dataAtPts->nbUniq;
837 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
838 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
840 double log_det = t_log_stretch_h1(
i,
i);
841 double log_det2 = log_det * log_det;
844 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
846 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
847 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
851 (A * 2 * nohat_k) * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j));
852 t_B_diff(
i,
j) = (
B * 2 * hat_k) * log_det *
t_kd(
i,
j) +
856 t_A_diff(
i,
j) * (t_dev(
m,
n) * t_diff_deviator(
m,
n,
k,
l))
860 A * t_diff_deviator(
m,
n,
i,
j) * t_diff_deviator(
m,
n,
k,
l)
866 t_dP(L,
J) = -t_L(
i,
j, L) *
873 (alphaU * ts_a) * (t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)
883 t_sym(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
884 t_approx_P_adjoint__dstretch(
j,
i));
889 t_dP(L,
J) += t_L(
i,
j, L) *
890 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
896 ++t_approx_P_adjoint__dstretch;
902 return getFTensor2FromMat<size_symm, size_symm>(dP);
905 auto t_dP = get_dP();
906 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
910 for (; rr != row_nb_dofs / 6; ++rr) {
912 auto t_m = get_ftensor2(K, 6 * rr, 0);
913 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
914 const double b =
a * t_row_base_fun * t_col_base_fun;
915 t_m(L,
J) -= b * t_dP(L,
J);
922 for (; rr != row_nb_base_functions; ++rr) {
933 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
934 boost::shared_ptr<double> total_energy_ptr)
936 totalEnergyPtr(total_energy_ptr) {
940 "dataAtPts is not allocated. Please set it before "
941 "using this operator.");
955 int nb_integration_pts = getGaussPts().size2();
957 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
960 auto &mat_d = dataAtPts->matD;
963 "wrong matD size, should be 6 by 6 but is %d by %d",
size_symm,
967 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(mat_d.data().data());
969 dataAtPts->energyAtPts.resize(nb_integration_pts,
false);
970 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
972 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
974 t_energy = 0.5 * (t_log_u(
i,
j) * (t_D(
i,
j,
k,
l) * t_log_u(
k,
l)));
980 if (totalEnergyPtr) {
981 auto t_w = getFTensor0IntegrationWeight();
982 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
983 double loc_energy = 0;
984 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
985 loc_energy += t_energy * t_w;
989 *totalEnergyPtr += getMeasure() * loc_energy;
996 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
997 boost::shared_ptr<MatrixDouble> strain_ptr,
998 boost::shared_ptr<MatrixDouble> stress_ptr,
999 boost::shared_ptr<HMHHencky> hencky_ptr)
1001 strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {}
1015 auto nb_integration_pts = stressPtr->size2();
1017 if (nb_integration_pts != getGaussPts().size2()) {
1019 "inconsistent number of integration points");
1023 auto get_D = [&]() {
1025 for (
auto &b : henckyPtr->blockData) {
1027 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
1028 dataAtPts->mu = b.shearModulusG;
1029 dataAtPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
1030 CHKERR henckyPtr->getMatDPtr(
1031 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
1032 dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
1033 CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(),
1034 b.bulkModulusK, b.shearModulusG);
1039 const auto E = henckyPtr->E;
1040 const auto nu = henckyPtr->nu;
1048 CHKERR henckyPtr->getMatDPtr(
1049 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
1058 strainPtr->resize(
size_symm, nb_integration_pts,
false);
1059 auto t_strain = getFTensor2SymmetricFromMat<3>(*strainPtr);
1060 auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressPtr);
1061 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
1062 &*dataAtPts->matInvD.data().begin());
1064 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
1065 &*dataAtPts->matD.data().begin());
1068 const double lambda = dataAtPts->lambda;
1069 const double mu = dataAtPts->mu;
1074 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1075 t_strain(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_stress(
k,
l);
1079 t_stress_symm_debug(
i,
j) = (t_stress(
i,
j) || t_stress(
j,
i)) / 2;
1081 t_stress_symm_debug_diff(
i,
j) =
1082 t_D(
i,
j,
k,
l) * t_strain(
k,
l) - t_stress_symm_debug(
i,
j);
1084 t_stress_symm_debug_diff(
i,
j) * t_stress_symm_debug_diff(
i,
j);
1085 double nrm0 = t_stress_symm_debug(
i,
j) * t_stress_symm_debug(
i,
j) +
1086 std::numeric_limits<double>::epsilon();
1087 constexpr double eps = 1e-10;
1088 if (std::fabs(std::sqrt(nrm / nrm0)) >
eps) {
1090 <<
"Stress symmetry check failed: " << std::endl
1091 << t_stress_symm_debug_diff << std::endl
1094 "Norm is too big: " + std::to_string(nrm / nrm0));
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
static PetscErrorCode ierr
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 ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
static constexpr auto size_symm
double young_modulus
Young modulus.
double poisson_ratio
Poisson ratio.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static double dynamicTime
static enum StretchSelector stretchSelector
static enum RotSelector gradApproximator
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
static PetscBool dynamicRelaxation
Calculate energy density for Hencky material model.
boost::shared_ptr< double > totalEnergyPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > strainPtr
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< HMHHencky > henckyPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< MatrixDouble > stressPtr
OpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
MoFEMErrorCode evaluateRhs(EntData &data)
boost::shared_ptr< HMHHencky > henckyPtr
MoFEMErrorCode evaluateLhs(EntData &data)
OpHenckyJacobian(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
OpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > &external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
MoFEMErrorCode integrate(EntData &data)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrateHencky(EntData &data)
MoFEMErrorCode integrate(EntData &data)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
VolUserDataOperator * returnOpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
MoFEMErrorCode extractBlockData(Sev sev)
Sev getOptionsSeverityLevels
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
VolUserDataOperator * returnOpCalculateVarStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
VolUserDataOperator * returnOpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
MoFEMErrorCode getMatDPtr(boost::shared_ptr< MatrixDouble > mat_D_ptr, boost::shared_ptr< MatrixDouble > mat_axiator_D_ptr, boost::shared_ptr< MatrixDouble > mat_deviator_D_ptr, double bulk_modulus_K, double shear_modulus_G)
MoFEMErrorCode getOptions(boost::shared_ptr< DataAtIntegrationPts > data_ptr)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEM::Interface & mField
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)
MoFEMErrorCode getInvMatDPtr(boost::shared_ptr< MatrixDouble > mat_inv_D_ptr, double bulk_modulus_K, double shear_modulus_G)
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.