1656 {
1658 if (type != MBVERTEX)
1660
1661 auto create_tag = [this](const std::string tag_name, const int size) {
1662 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1665 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1666 def_VAL);
1668 };
1669
1670 Tag th_w = create_tag("SpatialDisplacement", 3);
1671 Tag th_omega = create_tag("Omega", 3);
1672 Tag th_approxP = create_tag("Piola1Stress", 9);
1673 Tag th_sigma = create_tag("CauchyStress", 9);
1674 Tag th_log_u = create_tag("LogSpatialStreach", 9);
1675 Tag th_u = create_tag("SpatialStreach", 9);
1676 Tag th_h = create_tag("h", 9);
1677 Tag th_X = create_tag("X", 3);
1678 Tag th_detF = create_tag("detF", 1);
1679 Tag th_angular_momentum = create_tag("AngularMomentum", 3);
1680
1681 Tag th_u_eig_vec = create_tag("SpatialStreachEigenVec", 9);
1682 Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
1683
1684 int nb_gauss_pts = data.
getN().size1();
1685 if (
mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
1687 "Nb. of integration points is not equal to number points on "
1688 "post-processing mesh");
1689 }
1690
1691 auto t_w = getFTensor1FromMat<3>(
dataAtPts->wAtPts);
1692 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
1693 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
1694 auto t_log_u =
1695 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachTensorAtPts);
1696 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1697 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1698 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1699 auto t_coords = getFTensor1CoordsAtGaussPts();
1700
1705
1706
1709 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1713 &*
v.data().begin());
1715 };
1716
1719 &
m(0, 0), &
m(0, 1), &
m(0, 2),
1720
1721 &
m(1, 0), &
m(1, 1), &
m(1, 2),
1722
1723 &
m(2, 0), &
m(2, 1), &
m(2, 2));
1724
1725 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1727 t_m(
i,
j) = t_d(
i,
j);
1729 &*
m.data().begin());
1731 };
1732
1733 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1734
1735
1736 CHKERR save_vec_tag(th_w, t_w, gg);
1737 CHKERR save_vec_tag(th_X, t_coords, gg);
1738 CHKERR save_vec_tag(th_omega, t_omega, gg);
1739
1740
1741 CHKERR save_mat_tag(th_h, t_h, gg);
1742
1744 for (int ii = 0; ii != 3; ++ii)
1745 for (int jj = 0; jj != 3; ++jj)
1746 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
1747
1748 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
1749
1751 for (int ii = 0; ii != 3; ++ii)
1752 for (int jj = 0; jj != 3; ++jj)
1753 t_u_tmp(ii, jj) = t_u(ii, jj);
1754
1755 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
1756 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
1757
1760 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
1761 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
1763
1765 t_PhT(
i,
k) = t_approx_P(
i,
j) * t_R(
k,
j);
1768
1770 &t_leviPRT(0));
1771
1772 auto get_eiegn_vector_symmetric = [&](auto &t_u) {
1774
1775 for (int ii = 0; ii != 3; ++ii)
1776 for (int jj = 0; jj != 3; ++jj)
1777 t_m(ii, jj) = t_u(ii, jj);
1778
1780
1781
1782
1783 int n = 3, lda = 3, info, lwork = -1;
1784 double wkopt;
1786 &(eigen_values.data()[0]), &wkopt, lwork);
1787 if (info != 0)
1789 "is something wrong with lapack_dsyev info = %d", info);
1791 double work[lwork];
1793 &(eigen_values.data()[0]), work, lwork);
1794 if (info != 0)
1796 "is something wrong with lapack_dsyev info = %d", info);
1797
1799 &*
m.data().begin());
1801 &*eigen_values.data().begin());
1802
1804 };
1805
1806 CHKERR get_eiegn_vector_symmetric(t_u);
1807
1808 ++t_w;
1809 ++t_h;
1810 ++t_log_u;
1811 ++t_u;
1812 ++t_omega;
1813 ++t_R;
1814 ++t_approx_P;
1815 ++t_coords;
1816 }
1817
1819}
#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 ...
@ 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< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
VectorBoundedArray< double, 3 > VectorDouble3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
