1722 {
1724
1727 "Invalid data passed to Storakers stretch-from-stress operator");
1728 }
1730
1731 const int nb_integration_pts =
stressPtr->size1();
1732#ifndef NDEBUG
1733 if (nb_integration_pts != getGaussPts().size2()) {
1735 "Inconsistent number of Storakers integration points");
1736 }
1737#endif
1741 "Inconsistent number of external-pressure integration points");
1742 }
1743
1748 nb_integration_pts);
1750 DL>::size(*
dataAtPts->getDiffStretchTensorAtPts(),
1751 nb_integration_pts);
1753 *
dataAtPts->getEigenVals(), nb_integration_pts);
1755 DL>::size(*
dataAtPts->getEigenVecs(), nb_integration_pts);
1756 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
1760 DL>::size(
dataAtPts->matInvD, nb_integration_pts);
1762 DL>::size(
dataAtPts->logStretch2H1AtPts, nb_integration_pts);
1765 nb_integration_pts);
1767 DL>::size(
dataAtPts->rotMatAtPts, nb_integration_pts);
1770 nb_integration_pts);
1771
1773 dataAtPts->getStretchTensorAtPts()->clear();
1774 dataAtPts->getDiffStretchTensorAtPts()->clear();
1781 dataAtPts->logStretchTotalTensorAtPts.clear();
1782
1783 auto t_strain = getFTensor2SymmetricFromMat<
SPACE_DIM, -1,
DL>(*strainPtr);
1784 auto t_biot_stretch =
dataAtPts->getFTensorStretch(nb_integration_pts);
1785 auto t_diff_stretch =
dataAtPts->getFTensorDiffStretch(nb_integration_pts);
1787 auto t_omega =
dataAtPts->getFTensorRotAxis(nb_integration_pts);
1788 auto t_R =
dataAtPts->getFTensorRotMat(nb_integration_pts);
1789 auto t_biot_stress =
1790 dataAtPts->getFTensorAdjointPdstretch(nb_integration_pts);
1791 auto t_eigen_vals =
dataAtPts->getFTensorEigenVals(nb_integration_pts);
1792 auto t_eigen_vecs =
dataAtPts->getFTensorEigenVecs(nb_integration_pts);
1793 auto t_mat_d =
1795 auto t_mat_inv_d =
1797 auto t_log_u2_h1 =
dataAtPts->getFTensorLogStretch2H1(nb_integration_pts);
1798 auto t_log_stretch_total =
1799 dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
1802
1803 const auto parameters =
1804 storakersPtr->getMaterialParameters(getFEEntityHandle());
1805 StressToStrech stress_to_stretch(parameters);
1807 constexpr auto identity = [](const double value) { return value; };
1808
1811 "Rotation selector not handled by Storakers");
1812 }
1813
1814 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1816
1818 t_rotated_stress(
i,
j) = t_R(
k,
i) * t_stress(
k,
j);
1819 t_biot_stress(
i,
j) = t_diff_sym(
i,
j,
k,
l) * t_rotated_stress(
k,
l);
1820
1821 const double initial_y = stress_to_stretch.getScalarEquationState().y;
1823 const auto calculate_error = stress_to_stretch.calculate(
1824 t_biot_stress, t_biot_stretch, t_mat_d, t_mat_inv_d, q);
1825 if (calculate_error) {
1826
1827
1828 std::ostringstream atom_reproducer;
1829 atom_reproducer
1830 << std::scientific
1831 << std::setprecision(std::numeric_limits<double>::max_digits10)
1832 << "./storakers_stress_to_stretch_atom"
1833 << " -storakers_atom_eta " << parameters.eta
1834 << " -storakers_atom_mu " << parameters.mu
1835 << " -storakers_atom_beta " << parameters.beta
1836 << " -storakers_atom_q " << q
1837 << " -storakers_atom_initial_y " << initial_y
1838 << " -storakers_atom_solver_tol "
1839 << stress_to_stretch.absoluteTolerance
1840 << " -storakers_atom_biot_stress " << t_biot_stress(0, 0) << ","
1841 << t_biot_stress(0, 1) << "," << t_biot_stress(0, 2) << ","
1842 << t_biot_stress(1, 1) << "," << t_biot_stress(1, 2) << ","
1843 << t_biot_stress(2, 2) << " -log_no_color";
1844
1845 PetscMPIInt rank = -1;
1846 (void)MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1847 std::ostringstream atom_reproducer_log;
1848 atom_reproducer_log << "[" << rank
1849 << "] <warning> atom_reproducer_command="
1850 << atom_reproducer.str() << '\n';
1851 const auto atom_reproducer_log_string = atom_reproducer_log.str();
1852 std::fwrite(atom_reproducer_log_string.data(), 1,
1853 atom_reproducer_log_string.size(), stderr);
1854 std::fflush(stderr);
1855 }
1857
1858 const auto &state = stress_to_stretch.getScalarEquationState();
1859 const auto &eigenvectors = stress_to_stretch.getEigenVectors();
1860 t_eigen_vals(
i) = state.tLogU(
i);
1861 t_eigen_vecs(
i,
j) = eigenvectors(
i,
j);
1862 for (
int aa = 0; aa !=
SPACE_DIM; ++aa) {
1863 if (!std::isfinite(t_eigen_vals(aa))) {
1865 "Non-finite Storakers strain eigenvalue");
1866 }
1867 }
1868
1869 auto t_log_stretch =
1871 t_strain(
i,
j) = t_log_stretch(
i,
j);
1872 t_log_stretch_total(
i,
j) = t_strain(
i,
j);
1873 t_log_u2_h1(
i,
j) = 0;
1874
1875 t_nb_uniq =
1879 CHKERR sortEigenVals<SPACE_DIM>(
1883 }
1885 t_eigen_vals, t_eigen_vecs,
1886 static_cast<double (*)(double)>(std::exp),
1887 static_cast<double (*)(double)>(std::exp), t_nb_uniq);
1888 t_diff_stretch(
i,
j,
k,
l) = t_diff_stretch_mat(
i,
j,
k,
l);
1889
1890 ++t_strain;
1891 ++t_biot_stretch;
1892 ++t_diff_stretch;
1893 ++t_stress;
1894 ++t_omega;
1895 ++t_R;
1896 ++t_biot_stress;
1897 ++t_eigen_vals;
1898 ++t_eigen_vecs;
1899 ++t_mat_d;
1900 ++t_mat_inv_d;
1901 ++t_log_u2_h1;
1902 ++t_log_stretch_total;
1903 ++t_nb_uniq;
1904 }
1905
1907}
#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
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
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.
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
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
static enum RotSelector rotSelector
static auto exp(A &&t_w_vee, B &&theta)
