990 {
992
995
996 const int nb_integration_pts =
stressPtr->size1();
997#ifndef NDEBUG
998 if (nb_integration_pts != getGaussPts().size2()) {
1000 "inconsistent number of integration points");
1001 }
1002#endif
1006 "Inconsistent number of external-pressure integration points");
1007 }
1008
1013 nb_integration_pts);
1015 DL>::size(*
dataAtPts->getDiffStretchTensorAtPts(),
1016 nb_integration_pts);
1018 *
dataAtPts->getEigenVals(), nb_integration_pts);
1020 DL>::size(*
dataAtPts->getEigenVecs(), nb_integration_pts);
1021 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
1025 DL>::size(
dataAtPts->matInvD, nb_integration_pts);
1027 DL>::size(
dataAtPts->logStretch2H1AtPts, nb_integration_pts);
1030 nb_integration_pts);
1032 DL>::size(
dataAtPts->rotMatAtPts, nb_integration_pts);
1035 nb_integration_pts);
1036
1038 dataAtPts->getStretchTensorAtPts()->clear();
1039 dataAtPts->getDiffStretchTensorAtPts()->clear();
1046 dataAtPts->logStretchTotalTensorAtPts.clear();
1047
1048 auto t_strain = getFTensor2SymmetricFromMat<
SPACE_DIM, -1,
DL>(*strainPtr);
1049 auto t_biot_stretch =
dataAtPts->getFTensorStretch(nb_integration_pts);
1050 auto t_diff_stretch =
dataAtPts->getFTensorDiffStretch(nb_integration_pts);
1052 auto t_omega =
dataAtPts->getFTensorRotAxis(nb_integration_pts);
1053 auto t_R =
dataAtPts->getFTensorRotMat(nb_integration_pts);
1054 auto t_biot_stress =
1055 dataAtPts->getFTensorAdjointPdstretch(nb_integration_pts);
1056 auto t_eigen_vals =
dataAtPts->getFTensorEigenVals(nb_integration_pts);
1057 auto t_eigen_vecs =
dataAtPts->getFTensorEigenVecs(nb_integration_pts);
1058 auto t_mat_d =
1060 auto t_mat_inv_d =
1062 auto t_log_u2_h1 =
dataAtPts->getFTensorLogStretch2H1(nb_integration_pts);
1063 auto t_log_stretch_total =
1064 dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
1067
1068 const auto [def_c10, def_K] =
1070 const double c10 = def_c10;
1071 const double bulk_modulus = def_K;
1073 "no-stretch integration point");
1074
1077 "Rotation selector not handled by Abaqus Neo-Hookean");
1078 }
1079
1081 stretch_from_stress.
c10 = c10;
1082 stretch_from_stress.K = bulk_modulus;
1083 stretch_from_stress.allowNonConverged = true;
1085
1086 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1088
1090 t_rotated_stress(
i,
j) = t_R(
k,
i) * t_stress(
k,
j);
1091 t_biot_stress(
i,
j) = t_diff_sym(
i,
j,
k,
l) * t_rotated_stress(
k,
l);
1092 stretch_from_stress.tBiotStress(
i,
j) =
1093 (t_biot_stress(
i,
j) || t_biot_stress(
j,
i)) / 2.;
1094 stretch_from_stress.tHenckyStreach(
i,
j) = t_strain(
i,
j);
1095 stretch_from_stress.externalPressure =
1097
1098 CHKERR stretch_from_stress.calculateStretch();
1099 CHKERR stretch_from_stress.calculateDHenckyStreachDBiotStress();
1100
1101 t_strain(
i,
j) = stretch_from_stress.tHenckyStreach(
i,
j);
1102 t_log_stretch_total(
i,
j) = stretch_from_stress.tHenckyStreach(
i,
j);
1103 t_log_u2_h1(
i,
j) = 0;
1104 t_biot_stretch(
i,
j) = stretch_from_stress.tBiotStreach(
i,
j);
1105
1106 t_eigen_vals(
i) = stretch_from_stress.tHenckyStreachEigenVals(
i);
1107 t_eigen_vecs(
i,
j) = stretch_from_stress.tBiotStressEigenVecs(
i,
j);
1108 t_nb_uniq =
1112 CHKERR sortEigenVals<SPACE_DIM>(
1116 }
1117
1119 t_eigen_vals, t_eigen_vecs, static_cast<double (*)(double)>(std::exp),
1120 static_cast<double (*)(double)>(std::exp), t_nb_uniq);
1121 t_diff_stretch(
i,
j,
k,
l) = t_diff_stretch_mat(
i,
j,
k,
l);
1122
1124 auto get_d_biot_stretch_d_biot_stress = MatrixSizeHelper<
1126 DL>::size(d_biot_stretch_d_biot_stress_mat, 1);
1127 auto t_d_biot_stretch_d_biot_stress = get_d_biot_stretch_d_biot_stress();
1128
1129 t_d_biot_stretch_d_biot_stress(
i,
j,
k,
l) =
1130 t_diff_stretch(
i,
j,
m,
n) *
1131 stretch_from_stress.tDHenckyStreachDBiotStress(
m,
n,
k,
l);
1132
1133
1134
1135
1137 auto t_packed_compliance = getFTensor2FromMat<size_symm, size_symm>(
1138 d_biot_stretch_d_biot_stress_mat);
1140 t_compliance_eig_vecs(L,
J) =
1141 0.5 * (t_packed_compliance(L,
J) + t_packed_compliance(
J, L));
1144 t_compliance_eig_vals);
1145 const double numerical_relative_floor =
1146 std::sqrt(std::numeric_limits<double>::epsilon());
1147 double compliance_spectral_radius = 0.;
1148 for (
int aa = 0; aa !=
size_symm; ++aa) {
1149 compliance_spectral_radius = std::max(
1150 compliance_spectral_radius, std::abs(t_compliance_eig_vals(aa)));
1151 }
1152
1153
1154
1155
1156 const double material_compliance_eigenvalue_floor =
1157 stretch_from_stress.minimEigenValue /
1158 std::max(std::abs(stretch_from_stress.K),
1159 std::numeric_limits<double>::min());
1160 const double numerical_compliance_eigenvalue_floor =
1161 numerical_relative_floor * compliance_spectral_radius;
1162 const double compliance_eigenvalue_floor =
1163 std::max(material_compliance_eigenvalue_floor,
1164 numerical_compliance_eigenvalue_floor);
1165 bool project_compliance = false;
1166 for (
int aa = 0; aa !=
size_symm; ++aa) {
1167 project_compliance =
1168 project_compliance ||
1169 t_compliance_eig_vals(aa) < compliance_eigenvalue_floor;
1170 }
1171 if (project_compliance) {
1172 auto floor_compliance_eigenvalue =
1173 [compliance_eigenvalue_floor](
const double v) {
1174 return (
v + compliance_eigenvalue_floor +
1175 std::abs(
v - compliance_eigenvalue_floor)) /
1176 2.;
1177 };
1178 auto t_projected_compliance =
1180 floor_compliance_eigenvalue);
1181 t_packed_compliance(L,
J) = t_projected_compliance(L,
J);
1182 }
1183
1184
1185 auto d_biot_stress_d_hencky_mat = d_biot_stretch_d_biot_stress_mat;
1187 t_mat_inv_d(
i,
j,
k,
l) = t_d_biot_stretch_d_biot_stress(
i,
j,
k,
l);
1188
1189
1190
1191 for (
int row = 0; row !=
size_symm; ++row) {
1192 for (const int col : {1, 2, 4}) {
1193 d_biot_stress_d_hencky_mat(row, col) *= 2.;
1194 }
1195 }
1197 for (
int row = 0; row !=
size_symm; ++row) {
1198 for (const int col : {1, 2, 4}) {
1199 d_biot_stress_d_hencky_mat(row, col) *= 0.5;
1200 }
1201 }
1202 auto get_d_biot_stress_d_hencky = MatrixSizeHelper<
1204 DL>::size(d_biot_stress_d_hencky_mat, 1);
1205 auto t_d_biot_stress_d_hencky = get_d_biot_stress_d_hencky();
1206 t_mat_d(
i,
j,
k,
l) = t_d_biot_stress_d_hencky(
i,
j,
k,
l);
1207
1208 ++t_strain;
1209 ++t_biot_stretch;
1210 ++t_diff_stretch;
1211 ++t_stress;
1212 ++t_omega;
1213 ++t_R;
1214 ++t_biot_stress;
1215 ++t_eigen_vals;
1216 ++t_eigen_vecs;
1217 ++t_mat_d;
1218 ++t_mat_inv_d;
1219 ++t_log_u2_h1;
1220 ++t_log_stretch_total;
1221 ++t_nb_uniq;
1222 }
1223
1225}
#define FTENSOR_INDEXES(DIM,...)
Fourth-order symmetrization tensor.
#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< 'i', SPACE_DIM > i
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 getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
static constexpr auto size_symm
UBlasMatrix< double > MatrixDouble
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
FTensor::Index< 'm', 3 > m
static enum RotSelector rotSelector
static MoFEMErrorCode validateMaterialParameters(const double c10, const double K, const char *source)
static auto exp(A &&t_w_vee, B &&theta)
