diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 507ae8bdb..43e367ff8 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -16,11 +16,13 @@ namespace consts { -const int maxNSD = 3; +constexpr int maxNSD = 3; -const int maxNProp = 20; +constexpr int maxNoN = 27; // Max node count in nn_elem_props.h -const int maxOutput = 5; +constexpr int maxNProp = 20; + +constexpr int maxOutput = 5; /// Use inf numeric values to represent a value that is not set. const int int_inf = std::numeric_limits::infinity(); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index c0ef6eb65..9f151368c 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -1456,9 +1456,6 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double start_time = utils::cput(); #endif - // Maximum size of arrays sized by (3,eNoNw) -> (3,MAX_SIZE). - const int MAX_SIZE = 27; - using namespace consts; int cEq = com_mod.cEq; @@ -1583,7 +1580,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e es[1][2] = es[2][1]; es[2][0] = es[0][2]; - double esNx[3][MAX_SIZE]; + double esNx[3][maxNoN]; for (int a = 0; a < eNoNw; a++) { esNx[0][a] = es[0][0]*Nwx(0,a) + es[1][0]*Nwx(1,a) + es[2][0]*Nwx(2,a); @@ -1647,7 +1644,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Stabilization parameters // double up[3] = {}; - double updu[3][3][MAX_SIZE] = {}; + double updu[3][3][maxNoN] = {}; double tauM = 0.0; if (vmsFlag) { @@ -1782,9 +1779,6 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double start_time = utils::cput(); #endif - // Maximum size of arrays sized by (3,eNoNw) -> (3,MAX_SIZE). - const int MAX_SIZE = 27; - using namespace consts; int cEq = com_mod.cEq; @@ -1927,7 +1921,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e es[1][2] = es[2][1]; es[2][0] = es[0][2]; - double esNx[3][MAX_SIZE]; + double esNx[3][maxNoN]; for (int a = 0; a < eNoNw; a++) { esNx[0][a] = es[0][0]*Nwx(0,a) + es[1][0]*Nwx(1,a) + es[2][0]*Nwx(2,a); @@ -2100,10 +2094,10 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Local residual // - double updu[3][3][MAX_SIZE] = {}; - double uNx[MAX_SIZE] = {}; - double upNx[MAX_SIZE] = {}; - double uaNx[MAX_SIZE] = {}; + double updu[3][3][maxNoN] = {}; + double uNx[maxNoN] = {}; + double upNx[maxNoN] = {}; + double uaNx[maxNoN] = {}; for (int a = 0; a < eNoNw; a++) { lR(0,a) = lR(0,a) + wr*Nw(a)*rV[0] + w*(Nwx(0,a)*rM[0][0] + Nwx(1,a)*rM[1][0] + Nwx(2,a)*rM[2][0]); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 4b7d32c13..bfd54c254 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -102,6 +102,14 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So ya_l_s = 0.0; ya_l_n = 0.0; + if (lM.fN.size() != 0) { + for (int iFn = 0; iFn < nFn; iFn++) { + for (int i = 0; i < nsd; i++) { + fN(i,iFn) = lM.fN(i+nsd*iFn,e); + } + } + } + for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); ptr(a) = Ac; @@ -116,14 +124,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So dl(i,a) = Dg(i,Ac); } - if (lM.fN.size() != 0) { - for (int iFn = 0; iFn < nFn; iFn++) { - for (int i = 0; i < nsd; i++) { - fN(i,iFn) = lM.fN(i+nsd*iFn,e); - } - } - } - if (pS0.size() != 0) { pS0l.set_col(a, pS0.col(Ac)); } @@ -186,7 +186,11 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } } - if (g == 0 || !fs_1[0].lShpF) { + // Shape function gradients and the viscous response are constant + // within linear triangles and tetrahedra. + const bool recompute_visc = (g == 0 || !fs_1[0].lShpF); + + if (recompute_visc) { auto Nx = fs_1[0].Nx.rslice(g); nn::gnn(fs_1[0].eNoN, nsd, nsd, Nx, xwl, Nwx, Jac, ksix); if (utils::is_zero(Jac)) { @@ -220,7 +224,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N0 = fs_1[0].N.col(g); struct_ns::struct_3d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK); + ya_l_s, ya_l_n, lR, lK, recompute_visc); } break; case Equation_lElas: throw std::runtime_error("[construct_fsi] LELAS3D not implemented"); @@ -233,7 +237,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So ustruct::ustruct_3d_m(com_mod, cep_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, ya_l_n, lR, - lK, lKd); + lK, lKd, recompute_visc); break; } @@ -256,7 +260,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N0 = fs_1[0].N.col(g); struct_ns::struct_2d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK); + ya_l_s, ya_l_n, lR, lK, recompute_visc); } break; case Equation_ustruct: diff --git a/Code/Source/solver/mat_fun.h b/Code/Source/solver/mat_fun.h index f1b79b10d..c69782465 100644 --- a/Code/Source/solver/mat_fun.h +++ b/Code/Source/solver/mat_fun.h @@ -9,6 +9,7 @@ #include #include "Array.h" +#include "consts.h" #include "Tensor4.h" #include "Vector.h" #include "FE/Common/FEException.h" @@ -21,13 +22,79 @@ /// \todo [TODO:DaveP] this should just be a namespace? // namespace mat_fun { - // Define templated type aliases for Eigen matrices and tensors for convenience - template + /// @brief A 2nd order tensor, nsd x nsd, fixed size and stack allocated. + /// Used for the deformation gradient, stresses and similar quantities that + /// have a known size at compile time. + template using Matrix = Eigen::Matrix; - template + /// @brief A 4th order tensor, nsd x nsd x nsd x nsd, fixed size and stack + /// allocated. Used for the material elasticity tensor and other 4th order tensors + /// that have a known size at compile time. + template using Tensor = Eigen::TensorFixedSize>; + /// @brief One nsd-vector per element node, so nsd x eNoN. Row count is fixed + /// at compile time while column count is the element's node count, known only + /// at run time, so it is bounded by consts::maxNoN to stay stack allocated. + /// Used for shape function gradients and other per-node vector quantities. + template + using NodalMatrix = Eigen::Matrix; + + /// @brief One scalar per element node, so eNoN entries. Dynamic length bounded + /// by consts::maxNoN to stay stack allocated, as for NodalMatrix. Used for shape + /// function values and other per-node scalar quantities. + using NodalVector = Eigen::Matrix; + + // The eigen_view overloads below wrap an Array or Vector in an Eigen::Map that + // shares its storage, so the container must outlive the view. + + /// @brief Read-only Eigen view of an Array, sharing its storage. + /// + /// @tparam rows Row count, fixed at compile time; the columns are taken from the Array. + template + Eigen::Map> + eigen_view(const Array& A) { + if (A.nrows() != rows) { + svmp::raise( + "A view of " + std::to_string(rows) + " rows was requested for an array with " + + std::to_string(A.nrows()) + " rows."); + } + return {A.data(), rows, A.ncols()}; + } + + /// @brief Read-only Eigen view of a whole Array, sharing its storage. + inline Eigen::Map + eigen_view(const Array& A) { + return {A.data(), A.nrows(), A.ncols()}; + } + + /// @brief Writable Eigen view of a whole Array, sharing its storage. + inline Eigen::Map + eigen_view_mutable(Array& A) { + return {A.data(), A.nrows(), A.ncols()}; + } + + /// @brief Read-only Eigen view of a whole Vector, sharing its storage. + inline Eigen::Map + eigen_view(const Vector& v) { + return {v.data(), v.size()}; + } + + /// @brief Read-only Eigen view of a Vector, sharing its storage. + /// + /// @tparam rows Entry count, fixed at compile time. + template + Eigen::Map> + eigen_view(const Vector& v) { + if (v.size() != rows) { + svmp::raise( + "A view of " + std::to_string(rows) + " entries was requested for a vector with " + + std::to_string(v.size()) + " entries."); + } + return Eigen::Map>(v.data()); + } + // Function to convert Array to Eigen::Matrix template MatrixType convert_to_eigen_matrix(const Array& src) { @@ -179,14 +246,32 @@ namespace mat_fun { Tensor4 ten_ddot_3424(const Tensor4& A, const Tensor4& B, const int nd); /** - * @brief Contracts two 4th order tensors A and B over two dimensions, - * + * @brief Contracts two 4th order tensors A and B over two dimensions. + * + * For example, if dimsA = {0, 1} and dimsB = {2, 3} this is + * C_klmn = A_ijkl B_mnij (sum over i, j) + * + * @tparam nsd Number of spatial dimensions; each tensor is nsd^4. + * @param[in] A,B Fourth order tensors to contract. + * @param[in] dimsA,dimsB Indices of the contracted dimensions of A and B. + * @return The contracted tensor. */ template Tensor double_dot_product(const Tensor& A, const std::array& dimsA, const Tensor& B, const std::array& dimsB) { + // Fast path for dimsA = dimsB = {2,3}: C_ijmn = A_ijkl * B_mnkl. + if (dimsA[0] == 2 && dimsA[1] == 3 && dimsB[0] == 2 && dimsB[1] == 3) { + constexpr int N = nsd * nsd; + Tensor C; + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data()); + Eigen::Map> c(C.data()); + c.noalias() = a * b.transpose(); + return C; + } + // Define the contraction dimensions Eigen::array, 2> contractionDims = { Eigen::IndexPair(dimsA[0], dimsB[0]), // Contract A's dimsA[0] with B's dimsB[0] @@ -195,9 +280,6 @@ namespace mat_fun { // Return the double dot product return A.contract(B, contractionDims); - - // For some reason, in this case the Eigen::Tensor contract function is - // faster than a for loop implementation. } Tensor4 ten_dyad_prod(const Array& A, const Array& B, const int nd); @@ -215,20 +297,13 @@ namespace mat_fun { dyadic_product(const Matrix& A, const Matrix& B) { // Initialize the result tensor Tensor C; + constexpr int N = nsd * nsd; - // Compute the dyadic product: C_ijkl = A_ij * B_kl - for (int i = 0; i < nsd; ++i) { - for (int j = 0; j < nsd; ++j) { - for (int k = 0; k < nsd; ++k) { - for (int l = 0; l < nsd; ++l) { - C(i,j,k,l) = A(i,j) * B(k,l); - } - } - } - } - // For some reason, in this case the Eigen::Tensor contract function is - // slower than the for loop implementation - + // Column-major storage flattens index pairs: c(ij,kl) = a(ij) * b(kl). + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data()); + Eigen::Map> c(C.data()); + c.noalias() = a * b.transpose(); return C; } @@ -265,8 +340,9 @@ namespace mat_fun { /// @brief Create a 4th order tensor from symmetric outer product of two matrices: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk) /// - /// Reproduces 'FUNCTION TEN_SYMMPROD(A, B, nd) RESULT(C)'. - // + /// @tparam nsd Number of spatial dimensions. + /// @param[in] A,B Second order tensors. + /// @return The resulting 4th order tensor. template Tensor symmetric_dyadic_product(const Matrix& A, const Matrix& B) { @@ -275,17 +351,14 @@ namespace mat_fun { Tensor C; // Compute the symmetric product: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk) - for (int i = 0; i < nsd; ++i) { - for (int j = 0; j < nsd; ++j) { - for (int k = 0; k < nsd; ++k) { - for (int l = 0; l < nsd; ++l) { - C(i,j,k,l) = 0.5 * (A(i,k) * B(j,l) + A(i,l) * B(j,k)); - } - } + for (int l = 0; l < nsd; ++l) { + for (int k = 0; k < nsd; ++k) { + // blk views the (k,l) block of C, so blk(i,j) is C(i,j,k,l). + Eigen::Map> blk(C.data() + nsd * nsd * (k + nsd * l)); + blk.noalias() = 0.5 * (A.col(k) * B.col(l).transpose() + + A.col(l) * B.col(k).transpose()); } } - // For some reason, in this case the for loop implementation is faster - // than the Eigen::Tensor contract method // Return the symmetric product return C; diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index b23120176..7be520658 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -5,6 +5,7 @@ #include "mat_models.h" +#include "consts.h" #include "mat_fun.h" #include "utils.h" #include "ArtificialNeuralNetMaterial.h" @@ -14,13 +15,6 @@ namespace mat_models { -// Define templated type aliases for Eigen matrices and 4th order tensors for convenience -template -using Matrix = Eigen::Matrix; - -template -using Tensor = Eigen::TensorFixedSize>; - /// @brief Compute active component of deformation gradient tensor for @@ -228,7 +222,7 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC) * for Hyperelastic Isotropic and Anisotropic Materials" by Cheng and Zhang. * */ -template +template std::pair, Tensor> bar_to_iso( const Matrix& S_bar, const Tensor &CC_bar, const double J2d, const Matrix& C, const Matrix& Ci) @@ -240,7 +234,7 @@ std::pair, Tensor> bar_to_iso( double r1 = J2d * double_dot_product(C, S_bar) / nsd; // Compute isochoric 2nd Piola-Kirchhoff stress - auto S_iso = J2d*S_bar - r1*Ci; + const Matrix S_iso = J2d*S_bar - r1*Ci; // Compute isochoric material elasticity tensor Tensor PP = fourth_order_identity() - (1.0/nsd) * dyadic_product(Ci, C); // Important: using auto here causes tests to fail @@ -266,8 +260,8 @@ std::pair, Tensor> bar_to_iso( * @return Normalized sheet-normal direction vector. * @throws std::runtime_error if directions are parallel or if called in 2D. */ -template -Eigen::Matrix compute_sheet_normal(const Eigen::Matrix& fl) +template +Eigen::Matrix compute_sheet_normal(const FiberRef& fl) { using namespace mat_fun; @@ -287,10 +281,10 @@ Eigen::Matrix compute_sheet_normal(const Eigen::Matrix +template void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Matrix &F, const int nfd, - const Eigen::Matrix fl, + const FiberRef &fl, const double ya_f, const double ya_s, const double ya_n, Matrix &S, Matrix<3 * (nsd - 1)> &Dm, double &Ja) { using namespace consts; @@ -816,6 +810,15 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, cc_to_voigt_eigen(CC, Dm); } +// Explicitly instantiate compute_pk2cc for 2D and 3D. +template void compute_pk2cc<2>(const ComMod&, const CepMod&, const dmnType&, + const Matrix<2>&, const int, const FiberRef<2>&, + const double, const double, const double, Matrix<2>&, Matrix<3>&, double&); + +template void compute_pk2cc<3>(const ComMod&, const CepMod&, const dmnType&, + const Matrix<3>&, const int, const FiberRef<3>&, + const double, const double, const double, Matrix<3>&, Matrix<6>&, double&); + /** * @brief Get the 2nd Piola-Kirchhoff stress tensor and material elasticity tensor. * @@ -830,18 +833,13 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& if (nsd == 2) { // Copy deformation gradient to Eigen matrix - auto F_2D = mat_fun::convert_to_eigen_matrix(F); + const auto F_2D = mat_fun::convert_to_eigen_matrix>(F); - // Copy fiber directions to Eigen matrix - Eigen::Matrix fl_2D(2, nfd); - for (int i = 0; i < nfd; i++) { - fl_2D(0, i) = fl(0, i); - fl_2D(1, i) = fl(1, i); - } + const auto fl_2D = eigen_view<2>(fl); // Initialize stress and elasticity tensors - Eigen::Matrix2d S_2D = Eigen::Matrix2d::Zero(); - Eigen::Matrix3d Dm_2D = Eigen::Matrix3d::Zero(); + Matrix<2> S_2D = Matrix<2>::Zero(); + Matrix<3> Dm_2D = Matrix<3>::Zero(); // Call templated function compute_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl_2D, ya_f, ya_s, ya_n, S_2D, Dm_2D, Ja); @@ -852,20 +850,13 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& } else if (nsd == 3) { // Copy deformation gradient to Eigen matrix - auto F_3D = mat_fun::convert_to_eigen_matrix(F); - - // Copy fiber directions to Eigen matrix - Eigen::Matrix fl_3D(3, nfd); - for (int i = 0; i < nfd; i++) { - fl_3D(0, i) = fl(0, i); - fl_3D(1, i) = fl(1, i); - fl_3D(2, i) = fl(2, i); - } + const auto F_3D = mat_fun::convert_to_eigen_matrix>(F); + + const auto fl_3D = eigen_view<3>(fl); // Initialize stress and elasticity tensors - Eigen::Matrix3d S_3D = Eigen::Matrix3d::Zero(); - Eigen::Matrix Dm_3D; - Dm_3D.setZero(); + Matrix<3> S_3D = Matrix<3>::Zero(); + Matrix<6> Dm_3D = Matrix<6>::Zero(); // Call templated function compute_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl_3D, ya_f, ya_s, ya_n, S_3D, Dm_3D, Ja); @@ -1559,73 +1550,54 @@ void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, } } +namespace { + /** * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid * with a viscous pseudo-potential model. - * This is defined by a viscous pseuo-potential + * + * This is defined by a viscous pseudo-potential * Psi = mu/2 * tr(E_dot^2) * The viscous 2nd Piola-Kirchhoff stress is given by - * Svis = dPsi/dE_dot + * Svis = dPsi/dE_dot * = mu * E_dot * = mu * 1/2 * F^T * (grad(v) + grad(v)^T) * F * = mu * 1/2 * ( (F^T * Grad(v)) + (F^T * Grad(v))^T ) - * + * * @tparam nsd Number of spatial dimensions - * @param mu Solid viscosity parameter - * @param eNoN Number of nodes in an element - * @param Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) - * @param vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) - * @param F Deformation gradient matrix - * @param Svis Viscous 2nd Piola-Kirchhoff stress matrix - * @param Kvis_u Viscous tangent matrix contribution due to displacement - * @param Kvis_v Visous tangent matrix contribution due to velocity + * @param[in] mu Solid viscosity parameter + * @param[in] eNoN Number of nodes in an element + * @param[in] Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param[in] vx Velocity gradient matrix w.r.t. reference configuration coordinates (dv/dX) + * @param[in] F Deformation gradient matrix + * @param[out] Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param[out] Kvis_u Viscous tangent matrix contribution due to displacement + * @param[out] Kvis_v Viscous tangent matrix contribution due to velocity */ -void compute_visc_stress_potential(const double mu, const int eNoN, const Array& Nx, const Array& vx, const Array& F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { - - using namespace consts; - using namespace mat_fun; - using namespace utils; - - // Number of spatial dimensions - int nsd = F.nrows(); - - // Initialize Svis, Kvis_u, Kvis_v to zero - Svis = 0.0; - Kvis_u = 0.0; - Kvis_v = 0.0; +template +void compute_visc_stress_potential(const double mu, const int eNoN, const Array& Nx, + const Matrix& vx, const Matrix& F, + Matrix& Svis, Array3& Kvis_u, Array3& Kvis_v) { // Required intermediate terms for stress and tangent - auto Ft = transpose(F); - auto F_Ft = mat_mul(F, Ft); - auto Ft_vx = mat_mul(Ft, vx); - auto vxt = transpose(vx); - auto F_vxt = mat_mul(F, vxt); - - //double F_Nx[nsd][eNoN] = {0}, vx_Nx[nsd][eNoN] = {0}; - Array F_Nx(nsd,eNoN), vx_Nx(nsd,eNoN); - - for (int a = 0; a < eNoN; ++a) { - for (int i = 0; i < nsd; ++i) { - for (int j = 0; j < nsd; ++j) { - F_Nx(i,a) += F(i,j) * Nx(j,a); - vx_Nx(i,a) += vx(i,j) * Nx(j,a); - } - } - } + const Matrix F_Ft = F * F.transpose(); + const Matrix Ft_vx = F.transpose() * vx; + const Matrix F_vxt = F * vx.transpose(); + + // F_Nx(i,a) = sum_j F(i,j) * Nx(j,a), and likewise for vx. + const auto Nxm = eigen_view(Nx); + const NodalMatrix F_Nx = F * Nxm; + const NodalMatrix vx_Nx = vx * Nxm; // 2nd Piola-Kirchhoff stress due to viscosity // Svis = mu * 1/2 * ( (F^T * dv/dX) + (F^T * dv/dX)^T ) - Svis = mu * mat_symm(Ft_vx, nsd); + Svis = 0.5 * mu * (Ft_vx + Ft_vx.transpose()); // Tangent matrix contributions due to viscosity for (int b = 0; b < eNoN; ++b) { for (int a = 0; a < eNoN; ++a) { - double Nx_Nx = 0.0; - for (int i = 0; i < nsd; ++i) { - Nx_Nx += Nx(i,a) * Nx(i,b); - } + const double Nx_Nx = Nxm.col(a).dot(Nxm.col(b)); for (int i = 0; i < nsd; ++i) { for (int j = 0; j < nsd; ++j) { @@ -1641,84 +1613,65 @@ void compute_visc_stress_potential(const double mu, const int eNoN, const Array< /** * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid * with a Newtonian fluid-like viscosity model. + * * The viscous deviatoric Cauchy stress is given by * sigma_vis_dev = 2 * mu * d_dev * where d_dev = 1/2 * (grad(v) + grad(v)^T) - 1/3 * (div(v)) * I * The viscous 2nd Piola-Kirchhoff stress is given by a pull-back operation * Svis = 2 * mu * J * F^-1 * d_dev * F^-T - * - * Note, there is likely an error/bug in the tangent contributions that leads to suboptimal nonlinear convergence - * + * + * Note, there is likely an error/bug in the tangent contributions + * that leads to suboptimal nonlinear convergence. + * * @tparam nsd Number of spatial dimensions - * @param mu Solid viscosity parameter - * @param eNoN Number of nodes in an element - * @param Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) - * @param vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) - * @param F Deformation gradient matrix - * @param Svis Viscous 2nd Piola-Kirchhoff stress matrix - * @param Kvis_u Viscous tangent matrix contribution due to displacement - * @param Kvis_v Visous tangent matrix contribution due to velocity + * @param[in] mu Solid viscosity parameter + * @param[in] eNoN Number of nodes in an element + * @param[in] Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param[in] vx Velocity gradient matrix w.r.t. reference configuration coordinates (dv/dX) + * @param[in] F Deformation gradient matrix + * @param[out] Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param[out] Kvis_u Viscous tangent matrix contribution due to displacement + * @param[out] Kvis_v Viscous tangent matrix contribution due to velocity */ -void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array& Nx, const Array& vx, const Array& F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { - using namespace consts; - using namespace mat_fun; - using namespace utils; - - // Number of spatial dimensions - int nsd = F.nrows(); - - // Initialize Svis, Kvis_u, Kvis_v to zero - Svis = 0.0; - Kvis_u = 0.0; - Kvis_v = 0.0; +template +void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array& Nx, + const Matrix& vx, const Matrix& F, + Matrix& Svis, Array3& Kvis_u, Array3& Kvis_v) { // Get identity matrix, Jacobian, and F^-1 - auto Idm = mat_id(nsd); - auto J = mat_det(F, nsd); - auto Fi = mat_inv(F, nsd); - - // Required intermediate terms for stress and tangent - // vx_Fi: Velocity gradient in current configuration - auto vx_Fi = mat_mul(vx, Fi); - auto vx_Fi_symm = mat_symm(vx_Fi, nsd); - // ddev: Deviatoric part of rate of strain tensor - auto ddev = mat_dev(vx_Fi_symm, nsd); - //double Nx_Fi[nsd][eNoN] = {0}, ddev_Nx_Fi[nsd][eNoN] = {0}, vx_Fi_Nx_Fi[nsd][eNoN] = {0}; - Array Nx_Fi(nsd,eNoN), ddev_Nx_Fi(nsd,eNoN), vx_Fi_Nx_Fi(nsd,eNoN); - for (int a = 0; a < eNoN; ++a) { - for (int i = 0; i < nsd; ++i) { - for (int j = 0; j < nsd; ++j) { - Nx_Fi(i,a) += Nx(j,a) * Fi(j,i); - } - } - } - - mat_mul(ddev, Nx_Fi, ddev_Nx_Fi); - mat_mul(vx_Fi, Nx_Fi, vx_Fi_Nx_Fi); + const auto Idm = Matrix::Identity(); + const double J = F.determinant(); + const Matrix Fi = F.inverse(); + + // vx_Fi: Velocity gradient in current configuration + const Matrix vx_Fi = vx * Fi; + // d: rate of deformation tensor, the symmetric velocity gradient + const Matrix d = 0.5 * (vx_Fi + vx_Fi.transpose()); + // ddev: its deviatoric part + const Matrix ddev = d - (d.trace() / nsd) * Idm; + + // Nx_Fi(i,a) = sum_j Nx(j,a) * Fi(j,i), which is Fi^T * Nx. + const auto Nxm = eigen_view(Nx); + const NodalMatrix Nx_Fi = Fi.transpose() * Nxm; + const NodalMatrix ddev_Nx_Fi = ddev * Nx_Fi; + const NodalMatrix vx_Fi_Nx_Fi = vx_Fi * Nx_Fi; // 2nd Piola-Kirchhoff stress due to viscosity // Svis = 2 * mu * J * F^-1 * d_dev * F^-T - auto Fit = transpose(Fi); - auto ddev_Fit = mat_mul(ddev, Fit); - auto Fi_ddev_Fit = mat_mul(Fi, ddev_Fit); - Svis = 2.0 * mu * J * Fi_ddev_Fit; + Svis.noalias() = (2.0 * mu * J) * (Fi * ddev * Fi.transpose()); // Tangent matrix contributions due to viscosity - double r2d = 2.0 / nsd; + constexpr double r2d = 2.0 / nsd; for (int b = 0; b < eNoN; ++b) { for (int a = 0; a < eNoN; ++a) { - double Nx_Fi_Nx_Fi = 0.0; - for (int i = 0; i < nsd; ++i) { - Nx_Fi_Nx_Fi += Nx_Fi(i,a) * Nx_Fi(i,b); - } + const double Nx_Fi_Nx_Fi = Nx_Fi.col(a).dot(Nx_Fi.col(b)); for (int i = 0; i < nsd; ++i) { for (int j = 0; j < nsd; ++j) { int ii = i * nsd + j; // Derivative of the residual w.r.t displacement - Kvis_u(ii,a,b) = mu * J * (2.0 * + Kvis_u(ii,a,b) = mu * J * (2.0 * (ddev_Nx_Fi(i,a) * Nx_Fi(j,b) - ddev_Nx_Fi(i,b) * Nx_Fi(j,a)) - (Nx_Fi_Nx_Fi * vx_Fi(i,j) + Nx_Fi(i,b) * vx_Fi_Nx_Fi(j,a) - r2d * Nx_Fi(i,a) * vx_Fi_Nx_Fi(j,b))); @@ -1732,34 +1685,46 @@ void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array< } } +} // namespace -/** - * @brief Get the solid viscous PK2 stress and corresponding tangent matrix contributions - * Calls the appropriate function based on the viscosity type, either viscous - * pseudo-potential or Newtonian viscosity model. - * - * @tparam nsd Number of spatial dimensions - * @param[in] lDmn Domain object - * @param[in] eNoN Number of nodes in an element - * @param[in] Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) - * @param[in] vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) - * @param[in] F Deformation gradient matrix - * @param[out] Svis Viscous 2nd Piola-Kirchhoff stress matrix - * @param[out] Kvis_u Viscous tangent matrix contribution due to displacement - * @param[out] Kvis_v Viscous tangent matrix contribution due to velocity - */ -void compute_visc_stress_and_tangent(const dmnType& lDmn, const int eNoN, const Array& Nx, const Array& vx, const Array& F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { +/// @brief Dispatches to the viscous pseudo-potential or Newtonian model, or +/// zeroes the contributions when the domain has no viscosity model. +template +void ViscousResponse::update(const dmnType& lDmn, const int eNoN, + const Array& Nx, const Matrix& vx, const Matrix& F, + const bool recompute) { + + // Reuse current stored values. + if (!recompute) { + return; + } + + // The buffers only need resizing when the element node count changes. + if (Kvis_u_.ncols() != eNoN) { + Kvis_u_.resize(nsd*nsd, eNoN, eNoN); + Kvis_v_.resize(nsd*nsd, eNoN, eNoN); + } switch (lDmn.solid_visc.viscType) { case consts::SolidViscosityModelType::viscType_Newtonian: - compute_visc_stress_newtonian(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + compute_visc_stress_newtonian(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis_, Kvis_u_, Kvis_v_); break; case consts::SolidViscosityModelType::viscType_Potential: - compute_visc_stress_potential(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + compute_visc_stress_potential(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis_, Kvis_u_, Kvis_v_); + break; + + default: + // No viscosity model for this domain. + Svis_.setZero(); + Kvis_u_ = 0.0; + Kvis_v_ = 0.0; break; } } +// Instantiate the dimensions the solver supports. +template class ViscousResponse<2>; +template class ViscousResponse<3>; + }; diff --git a/Code/Source/solver/mat_models.h b/Code/Source/solver/mat_models.h index 6a959e132..a7d0e2892 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -17,19 +17,27 @@ namespace mat_models { -void actv_strain(const ComMod& com_mod, const CepMod& cep_mod, const double gf, +using mat_fun::Matrix; +using mat_fun::Tensor; +using mat_fun::NodalMatrix; +using mat_fun::eigen_view; + +void actv_strain(const ComMod& com_mod, const CepMod& cep_mod, const double gf, const int nfd, const Array& fl, Array& Fa); void cc_to_voigt(const int nsd, const Tensor4& CC, Array& Dm); void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); +/// @brief Fiber directions at a Gauss point, one direction per column +/// with row count fixed at compile time to the number of spatial dimensions. +template +using FiberRef = Eigen::Ref>; + /** * @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors * including both dilational and isochoric components. * - * Reproduces the Fortran 'GETPK2CC' subroutine. - * * @param[in] com_mod Object containing global common variables. * @param[in] cep_mod Object containing electrophysiology-specific common * variables. @@ -45,7 +53,17 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); * @param[out] Ja Jacobian for active strain * * @return None, but modifies S, Dm, and Ja in place. + * + * @tparam nsd Number of spatial dimensions. */ +template +void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, + const dmnType &lDmn, const Matrix &F, const int nfd, + const FiberRef &fl, const double ya_f, + const double ya_s, const double ya_n, Matrix &S, + Matrix<3 * (nsd - 1)> &Dm, double &Ja); + +/// @brief Array-based overload, for callers whose dimension is a run-time value. void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Array &F, const int nfd, const Array &fl, const double ya_f, @@ -66,14 +84,47 @@ void compute_svol_p(const ComMod& com_mod, const CepMod& cep_mod, const stModelT void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, double& ro, double& bt, double& dro, double& dbt, const double Ja); -void compute_visc_stress_potential(const double mu, const int eNoN, const Array& Nx, const double vx, const double F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v); -void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array& Nx, const Array& vx, const Array& F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v); - -void compute_visc_stress_and_tangent(const dmnType& lDmn, const int eNoN, const Array& Nx, const Array& vx, const Array& F, - Array& Svis, Array3& Kvis_u, Array3& Kvis_v); +/// @brief Viscous 2nd Piola-Kirchhoff stress and tangent contributions at a +/// Gauss point, for the viscosity model configured for the domain. +/// +/// Owns the buffers it computes into, so a kernel can hold one instance across +/// an element's Gauss points and skip the calls it does not need. +/// +/// @tparam nsd Number of spatial dimensions. +/// +template +class ViscousResponse { + public: + /// @brief Evaluate the domain's viscosity model at this Gauss point. + /// + /// Does nothing when recompute is false, preserving values from the + /// previous call. Used for linear triangules and tetrahedra for which + /// Nx, vx, and F are constant across the element's Gauss points. + /// + /// @param[in] lDmn Domain, supplying the viscosity model and its parameters. + /// @param[in] eNoN Number of element nodes. + /// @param[in] Nx Shape function spatial derivatives. + /// @param[in] vx Velocity gradient. + /// @param[in] F Deformation gradient. + /// @param[in] recompute Whether the stored values must be recomputed. + void update(const dmnType& lDmn, const int eNoN, const Array& Nx, + const Matrix& vx, const Matrix& F, const bool recompute); + + /// @brief Viscous 2nd Piola-Kirchhoff stress. + const Matrix& S() const { return Svis_; } + + /// @brief Tangent w.r.t. displacement. du(i*nsd + j, a, b) is the (i,j) + /// entry of the block coupling nodes a and b. + double du(const int ij, const int a, const int b) const { return Kvis_u_(ij, a, b); } + + /// @brief Tangent w.r.t. velocity. + double dv(const int ij, const int a, const int b) const { return Kvis_v_(ij, a, b); } + + private: + Matrix Svis_; + Array3 Kvis_u_, Kvis_v_; +}; }; #endif diff --git a/Code/Source/solver/nn_elem_props.h b/Code/Source/solver/nn_elem_props.h index 7dc89e380..f63cc9243 100644 --- a/Code/Source/solver/nn_elem_props.h +++ b/Code/Source/solver/nn_elem_props.h @@ -25,7 +25,7 @@ SetElementPropsMapType set_3d_element_props = { mesh.nG = 6; mesh.vtkType = 13; mesh.nEf = 3; - mesh.lShpF = true; + mesh.lShpF = false; } }, diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index f350c8c6c..ac198c3ef 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -16,6 +16,7 @@ #include "nn.h" #include "utils.h" #include "DebugMsg.h" +#include namespace struct_ns { @@ -179,8 +180,12 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve } } -/// @brief Replicates the Fortan 'CONSTRUCT_dSOLID' subroutine. -// +/// @brief Assemble the residual and tangent contributions of one solid mesh. +/// +/// @param[in,out] com_mod Global common variables. +/// @param[in] cep_mod Electrophysiology variables, supplying the active stress. +/// @param[in] lM Mesh whose elements are assembled. +/// @param[in] solutions Acceleration, velocity and displacement. void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const SolutionStates& solutions) { const auto& Ag = solutions.intermediate.get_acceleration(); @@ -251,6 +256,14 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const ya_l_s = 0.0; ya_l_n = 0.0; + if (lM.fN.size() != 0) { + for (int iFn = 0; iFn < nFn; iFn++) { + for (int i = 0; i < nsd; i++) { + fN(i,iFn) = lM.fN(i+nsd*iFn,e); + } + } + } + for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); ptr(a) = Ac; @@ -266,14 +279,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const yl(i,a) = Yg(i,Ac); } - if (lM.fN.size() != 0) { - for (int iFn = 0; iFn < nFn; iFn++) { - for (int i = 0; i < nsd; i++) { - fN(i,iFn) = lM.fN(i+nsd*iFn,e); - } - } - } - if (pS0.size() != 0) { pS0l.set_col(a, pS0.col(Ac)); } @@ -294,7 +299,11 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array ksix(nsd,nsd); for (int g = 0; g < lM.nG; g++) { - if (g == 0 || !lM.lShpF) { + // Shape function gradients and the viscous response are constant + // within linear triangles and tetrahedra. + const bool recompute_visc = (g == 0 || !lM.lShpF); + + if (recompute_visc) { auto Nx_g = lM.Nx.slice(g); nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix); if (utils::is_zero(Jac)) { @@ -307,7 +316,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK); + pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK, recompute_visc); #if 0 if (e == 0 && g == 0) { @@ -321,16 +330,16 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } else if (nsd == 2) { struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK); + pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK, recompute_visc); } // Prestress if (pstEq) { for (int a = 0; a < eNoN; a++) { int Ac = ptr(a); - pSa(Ac) = pSa(Ac) + w*N(a); + pSa(Ac) += w*N(a); for (int i = 0; i < pSn.nrows(); i++) { - pSn(i,Ac) = pSn(i,Ac) + w*N(a)*pSl(i); + pSn(i,Ac) += w*N(a)*pSl(i); } } } @@ -349,7 +358,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &fN, const Array &pS0l, Vector &pSl, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK) { + Array &lR, Array3 &lK, const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -370,59 +379,45 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // double rho = dmn.prop.at(PhysicalPropertyType::solid_density); double dmp = dmn.prop.at(PhysicalPropertyType::damping); - Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), dmn.prop.at(PhysicalPropertyType::f_y)}); + const Eigen::Vector2d fb{dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y)}; double afu = eq.af * eq.beta*dt*dt; double afv = eq.af * eq.gam*dt; double amd = eq.am * rho + eq.af * eq.gam * dt * dmp; double afl = eq.af * eq.beta * dt * dt; int i = eq.s; - int j = i + 1; #ifdef debug_struct_2d dmsg << "i: " << i; - dmsg << "j: " << j; dmsg << "amd: " << amd; dmsg << "afl: " << afl; dmsg << "w: " << w; #endif - // Inertia, body force and deformation tensor (F) - // - Array F(2,2), S0(2,2), vx(2,2); - Vector ud(2); - - ud = -rho*fb; - F = 0.0; - F(0,0) = 1.0; - F(1,1) = 1.0; - S0 = 0.0; - - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - - for (int a = 0; a < eNoN; a++) { - ud(0) = ud(0) + N(a)*(rho*(al(i,a)-bfl(0,a)) + dmp*yl(i,a)); - ud(1) = ud(1) + N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a)); - - vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); - - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - - S0(0,0) = S0(0,0) + N(a)*pS0l(0,a); - S0(1,1) = S0(1,1) + N(a)*pS0l(1,a); - S0(0,1) = S0(0,1) + N(a)*pS0l(2,a); - - ya_g_f = ya_g_f + N(a) * ya_l_f(a); - ya_g_s = ya_g_s + N(a) * ya_l_s(a); - ya_g_n = ya_g_n + N(a) * ya_l_n(a); - } + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nxm = eigen_view<2>(Nx); // grad(N_a) per column + const auto Nm = eigen_view(N); // shape functions + const auto disp = eigen_view(dl).middleRows<2>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<2>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<2>(i); // nodal accelerations + const auto bfm = eigen_view<2>(bfl); // nodal body force + auto lRv = eigen_view_mutable(lR).topRows<2>(); // rows this kernel adds to + + // Inertia, damping and body force: the term the residual weights with N + const Eigen::Vector2d ud = (rho*(acc - bfm) + dmp*vel) * Nm - rho * fb; + + // Active stress activation along fiber, sheet and sheet-normal + const double ya_g_f = eigen_view(ya_l_f).dot(Nm); + const double ya_g_s = eigen_view(ya_l_s).dot(Nm); + const double ya_g_n = eigen_view(ya_l_n).dot(Nm); + + // Prestress at this Gauss point, in Voigt order [11, 22, 12] + const Eigen::Vector pS0g = eigen_view<3>(pS0l) * Nm; + + Matrix<2> S0; + S0 << pS0g(0), pS0g(2), + pS0g(2), pS0g(1); + #ifdef debug_struct_2d dmsg << "ud: " << ud(0) << " " << ud(1); dmsg << "F: " << F(0,0); @@ -431,23 +426,25 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, dmsg << "ya_g_n: " << ya_g_n; #endif - S0(1,0) = S0(0,1); + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<2> vx = vel * Nxm.transpose(); + const Matrix<2> F = Matrix<2>::Identity() + disp * Nxm.transpose(); // 2nd Piola-Kirchhoff stress (S) and material stiffness tensor in Voight notation (Dm) - Array S(2,2), Dm(3,3); + Matrix<2> S; + Matrix<3> Dm; double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, ya_g_s, + mat_models::compute_pk2cc<2>(com_mod, cep_mod, dmn, F, nFn, eigen_view<2>(fN), ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); - // Viscous 2nd Piola-Kirchhoff stress and tangent contributions - Array Svis(2,2); - Array3 Kvis_u(4, eNoN, eNoN); - Array3 Kvis_v(4, eNoN, eNoN); - - mat_models::compute_visc_stress_and_tangent(dmn, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. + // Reuse from the previous Gauss point when shape function gradients + // are constant within an element (e.g. linear triangles, tetrahedra). + static mat_models::ViscousResponse<2> visc; + visc.update(dmn, eNoN, Nx, vx, F, recompute_visc); // Elastic + Viscous stresses - S = S + Svis; + S = S + visc.S(); // Prestress pSl(0) = S(0,0); @@ -459,80 +456,60 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // 1st Piola-Kirchhoff tensor (P) // - Array P(2,2), DBm(3,2); - Array3 Bm(3,2,eNoN); - P = mat_fun::mat_mul(F, S); + const Matrix<2> P = F * S; #ifdef debug_struct_2d dmsg << "P: " << P(0,0) << " " << P(0,1); dmsg << " " << P(1,0) << " " << P(1,1); #endif - // Local residual - for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*ud(0) + Nx(0,a)*P(0,0) + Nx(1,a)*P(0,1)); - lR(1,a) = lR(1,a) + w*(N(a)*ud(1) + Nx(0,a)*P(1,0) + Nx(1,a)*P(1,1)); - } + // Local residual: inertia and body force, plus div P + lRv += w * (ud * Nm.transpose() + P * Nxm); - // Auxilary quantities for computing stiffness tensor + // Strain-displacement matrix; Bm[a] maps node a to Voigt strain // - for (int a = 0; a < eNoN; a++) { - Bm(0,0,a) = Nx(0,a)*F(0,0); - Bm(0,1,a) = Nx(0,a)*F(1,0); + std::array, consts::maxNoN> Bm; + const Matrix<2> Ft = F.transpose(); - Bm(1,0,a) = Nx(1,a)*F(0,1); - Bm(1,1,a) = Nx(1,a)*F(1,1); + for (int a = 0; a < eNoN; a++) { + const auto g = Nxm.col(a); // grad(N_a) - Bm(2,0,a) = (Nx(0,a)*F(0,1) + F(0,0)*Nx(1,a)); - Bm(2,1,a) = (Nx(0,a)*F(1,1) + F(1,0)*Nx(1,a)); + Bm[a].row(0) = g(0) * Ft.row(0); // dE_11 + Bm[a].row(1) = g(1) * Ft.row(1); // dE_22 + Bm[a].row(2) = g(0) * Ft.row(1) + g(1) * Ft.row(0); // 2 dE_12 } - Array NxFi(2,eNoN), DdNx(2,eNoN), VxNx(2,eNoN); - // Local stiffness tensor - double T1, NxNx, NxSNx, BmDBm; + double T1, NxSNx, BmDBm; - for (int b = 0; b < eNoN; b++) { - for (int a = 0; a < eNoN; a++) { - - // Geometric stiffness - NxSNx = Nx(0,a)*S(0,0)*Nx(0,b) + Nx(1,a)*S(1,0)*Nx(0,b) + - Nx(0,a)*S(0,1)*Nx(1,b) + Nx(1,a)*S(1,1)*Nx(1,b); - T1 = amd*N(a)*N(b) + afu*NxSNx; + for (int b = 0; b < eNoN; b++) { - // Material stiffness (Bt*D*B) - DBm(0,0) = Dm(0,0)*Bm(0,0,b) + Dm(0,1)*Bm(1,0,b) + Dm(0,2)*Bm(2,0,b); - DBm(0,1) = Dm(0,0)*Bm(0,1,b) + Dm(0,1)*Bm(1,1,b) + Dm(0,2)*Bm(2,1,b); + // Material stiffness for node b + const Eigen::Matrix DBm = Dm * Bm[b]; - DBm(1,0) = Dm(1,0)*Bm(0,0,b) + Dm(1,1)*Bm(1,0,b) + Dm(1,2)*Bm(2,0,b); - DBm(1,1) = Dm(1,0)*Bm(0,1,b) + Dm(1,1)*Bm(1,1,b) + Dm(1,2)*Bm(2,1,b); + // Geometric stiffness: S*grad(N_b) + const Eigen::Vector2d SNx = S * Nxm.col(b); - DBm(2,0) = Dm(2,0)*Bm(0,0,b) + Dm(2,1)*Bm(1,0,b) + Dm(2,2)*Bm(2,0,b); - DBm(2,1) = Dm(2,0)*Bm(0,1,b) + Dm(2,1)*Bm(1,1,b) + Dm(2,2)*Bm(2,1,b); + for (int a = 0; a < eNoN; a++) { + // Geometric stiffness + NxSNx = Nxm.col(a).dot(SNx); + T1 = amd*N(a)*N(b) + afu*NxSNx; // dM1/du1 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + Bm(2,0,a)*DBm(2,0); - - lK(0,a,b) = lK(0,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(0,a,b)) + afv*Kvis_v(0,a,b) ); + BmDBm = Bm[a].col(0).dot(DBm.col(0)); + lK(0,a,b) += w*( T1 + afu*(BmDBm + visc.du(0,a,b)) + afv*visc.dv(0,a,b) ); // dM1/du2 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,0,a)*DBm(0,1) + Bm(1,0,a)*DBm(1,1) + Bm(2,0,a)*DBm(2,1); - - lK(1,a,b) = lK(1,a,b) + w*( afu*(BmDBm + Kvis_u(1,a,b)) + afv*Kvis_v(1,a,b) ); + BmDBm = Bm[a].col(0).dot(DBm.col(1)); + lK(1,a,b) += w*( afu*(BmDBm + visc.du(1,a,b)) + afv*visc.dv(1,a,b) ); // dM2/du1 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,1,a)*DBm(0,0) + Bm(1,1,a)*DBm(1,0) + Bm(2,1,a)*DBm(2,0); - - lK(dof+0,a,b) = lK(dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(2,a,b)) + afv*Kvis_v(2,a,b) ); + BmDBm = Bm[a].col(1).dot(DBm.col(0)); + lK(dof+0,a,b) += w*( afu*(BmDBm + visc.du(2,a,b)) + afv*visc.dv(2,a,b) ); // dM2/du2 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1); - - lK(dof+1,a,b) = lK(dof+1,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(3,a,b)) + afv*Kvis_v(3,a,b) ); + BmDBm = Bm[a].col(1).dot(DBm.col(1)); + lK(dof+1,a,b) += w*( T1 + afu*(BmDBm + visc.du(3,a,b)) + afv*visc.dv(3,a,b) ); } } } @@ -545,7 +522,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &fN, const Array &pS0l, Vector &pSl, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK) { + Array &lR, Array3 &lK, const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -568,9 +545,9 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // double rho = dmn.prop.at(PhysicalPropertyType::solid_density); double dmp = dmn.prop.at(PhysicalPropertyType::damping); - Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), - dmn.prop.at(PhysicalPropertyType::f_y), - dmn.prop.at(PhysicalPropertyType::f_z)}); + const Eigen::Vector3d fb{dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y), + dmn.prop.at(PhysicalPropertyType::f_z)}; double afu = eq.af * eq.beta*dt*dt; double afv = eq.af * eq.gam*dt; @@ -585,88 +562,53 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, #endif int i = eq.s; - int j = i + 1; - int k = j + 1; - // Inertia, body force and deformation tensor (F) - // - Array F(3,3), S0(3,3), vx(3,3); - Vector ud(3); - - double F_f[3][3]={}; - F_f[0][0] = 1.0; - F_f[1][1] = 1.0; - F_f[2][2] = 1.0; + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nxm = eigen_view<3>(Nx); // grad(N_a) per column + const auto Nm = eigen_view(N); // shape functions + const auto disp = eigen_view(dl).middleRows<3>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<3>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<3>(i); // nodal accelerations + const auto bfm = eigen_view<3>(bfl); // nodal body force + auto lRv = eigen_view_mutable(lR).topRows<3>(); // rows this kernel adds to - ud = -rho*fb; - F = 0.0; - F(0,0) = 1.0; - F(1,1) = 1.0; - F(2,2) = 1.0; - S0 = 0.0; - - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - - for (int a = 0; a < eNoN; a++) { - ud(0) = ud(0) + N(a)*(rho*(al(i,a)-bfl(0,a)) + dmp*yl(i,a)); - ud(1) = ud(1) + N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a)); - ud(2) = ud(2) + N(a)*(rho*(al(k,a)-bfl(2,a)) + dmp*yl(k,a)); - - vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); - vx(0,2) = vx(0,2) + Nx(2,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); - vx(1,2) = vx(1,2) + Nx(2,a)*yl(j,a); - vx(2,0) = vx(2,0) + Nx(0,a)*yl(k,a); - vx(2,1) = vx(2,1) + Nx(1,a)*yl(k,a); - vx(2,2) = vx(2,2) + Nx(2,a)*yl(k,a); + // Inertia, damping and body force. + const Eigen::Vector3d ud = (rho*(acc - bfm) + dmp*vel) * Nm - rho * fb; - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); + // Active stress activation along fiber, sheet and sheet-normal + const double ya_g_f = eigen_view(ya_l_f).dot(Nm); + const double ya_g_s = eigen_view(ya_l_s).dot(Nm); + const double ya_g_n = eigen_view(ya_l_n).dot(Nm); - S0(0,0) = S0(0,0) + N(a)*pS0l(0,a); - S0(1,1) = S0(1,1) + N(a)*pS0l(1,a); - S0(2,2) = S0(2,2) + N(a)*pS0l(2,a); - S0(0,1) = S0(0,1) + N(a)*pS0l(3,a); - S0(1,2) = S0(1,2) + N(a)*pS0l(4,a); - S0(2,0) = S0(2,0) + N(a)*pS0l(5,a); + // Prestress at this Gauss point, in Voigt order [11, 22, 33, 12, 23, 31] + const Eigen::Vector pS0g = eigen_view<6>(pS0l) * Nm; - ya_g_f = ya_g_f + N(a) * ya_l_f(a); - ya_g_s = ya_g_s + N(a) * ya_l_s(a); - ya_g_n = ya_g_n + N(a) * ya_l_n(a); - } + Matrix<3> S0; + S0 << pS0g(0), pS0g(3), pS0g(5), + pS0g(3), pS0g(1), pS0g(4), + pS0g(5), pS0g(4), pS0g(2); - S0(1,0) = S0(0,1); - S0(2,1) = S0(1,2); - S0(0,2) = S0(2,0); + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<3> vx = vel * Nxm.transpose(); + const Matrix<3> F = Matrix<3>::Identity() + disp * Nxm.transpose(); // 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in - // Voigt notationa (Dm) + // Voigt notation (Dm) // - Array S(3,3), Dm(6,6); + Matrix<3> S; + Matrix<6> Dm; double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, ya_g_s, + mat_models::compute_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, eigen_view<3>(fN), ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); - // Viscous 2nd Piola-Kirchhoff stress and tangent contributions - Array Svis(3,3); - Array3 Kvis_u(9, eNoN, eNoN); - Array3 Kvis_v(9, eNoN, eNoN); - - mat_models::compute_visc_stress_and_tangent(dmn, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. + // Reuse from the previous Gauss point when shape function gradients + // are constant within an element (e.g. linear triangles, tetrahedra). + static mat_models::ViscousResponse<3> visc; + visc.update(dmn, eNoN, Nx, vx, F, recompute_visc); // Elastic + Viscous stresses - S = S + Svis; + S = S + visc.S(); #ifdef debug_struct_3d dmsg << "Jac: " << Jac; @@ -689,141 +631,78 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // 1st Piola-Kirchhoff tensor (P) // - Array P(3,3); - Array3 Bm(6,3,eNoN); - mat_fun::mat_mul(F, S, P); + const Matrix<3> P = F * S; - // Local residual - for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*ud(0) + Nx(0,a)*P(0,0) + Nx(1,a)*P(0,1) + Nx(2,a)*P(0,2)); - lR(1,a) = lR(1,a) + w*(N(a)*ud(1) + Nx(0,a)*P(1,0) + Nx(1,a)*P(1,1) + Nx(2,a)*P(1,2)); - lR(2,a) = lR(2,a) + w*(N(a)*ud(2) + Nx(0,a)*P(2,0) + Nx(1,a)*P(2,1) + Nx(2,a)*P(2,2)); - } + // Local residual: inertia and body force, plus div P + lRv += w * (ud * Nm.transpose() + P * Nxm); - // Auxilary quantities for computing stiffness tensor + // Strain-displacement matrix; Bm[a] maps node a to Voigt strain // - for (int a = 0; a < eNoN; a++) { - Bm(0,0,a) = Nx(0,a)*F(0,0); - Bm(0,1,a) = Nx(0,a)*F(1,0); - Bm(0,2,a) = Nx(0,a)*F(2,0); - - Bm(1,0,a) = Nx(1,a)*F(0,1); - Bm(1,1,a) = Nx(1,a)*F(1,1); - Bm(1,2,a) = Nx(1,a)*F(2,1); - - Bm(2,0,a) = Nx(2,a)*F(0,2); - Bm(2,1,a) = Nx(2,a)*F(1,2); - Bm(2,2,a) = Nx(2,a)*F(2,2); - - Bm(3,0,a) = (Nx(0,a)*F(0,1) + F(0,0)*Nx(1,a)); - Bm(3,1,a) = (Nx(0,a)*F(1,1) + F(1,0)*Nx(1,a)); - Bm(3,2,a) = (Nx(0,a)*F(2,1) + F(2,0)*Nx(1,a)); + std::array, consts::maxNoN> Bm; + const Matrix<3> Ft = F.transpose(); - Bm(4,0,a) = (Nx(1,a)*F(0,2) + F(0,1)*Nx(2,a)); - Bm(4,1,a) = (Nx(1,a)*F(1,2) + F(1,1)*Nx(2,a)); - Bm(4,2,a) = (Nx(1,a)*F(2,2) + F(2,1)*Nx(2,a)); - - Bm(5,0,a) = (Nx(2,a)*F(0,0) + F(0,2)*Nx(0,a)); - Bm(5,1,a) = (Nx(2,a)*F(1,0) + F(1,2)*Nx(0,a)); - Bm(5,2,a) = (Nx(2,a)*F(2,0) + F(2,2)*Nx(0,a)); + for (int a = 0; a < eNoN; a++) { + const auto g = Nxm.col(a); // grad(N_a) + + Bm[a].row(0) = g(0) * Ft.row(0); // dE_11 + Bm[a].row(1) = g(1) * Ft.row(1); // dE_22 + Bm[a].row(2) = g(2) * Ft.row(2); // dE_33 + Bm[a].row(3) = g(0) * Ft.row(1) + g(1) * Ft.row(0); // 2 dE_12 + Bm[a].row(4) = g(1) * Ft.row(2) + g(2) * Ft.row(1); // 2 dE_23 + Bm[a].row(5) = g(2) * Ft.row(0) + g(0) * Ft.row(2); // 2 dE_31 } // Local stiffness tensor - double NxSNx, T1, NxNx, BmDBm, Tv; - - Array DBm(6,3); + double NxSNx, T1, BmDBm; for (int b = 0; b < eNoN; b++) { - // Material stiffness (D*B). Shapes are fixed by the declarations above -- - // Dm(6,6), Bm(6,3,eNoN), DBm(6,3) -- so state them and skip the run-time - // shape check that the unparameterised overload would otherwise repeat on - // every one of these calls. - mat_mul<6, 6, 3>(Dm, Bm.rslice(b), DBm); + // Material stiffness for node b + const Eigen::Matrix DBm = Dm * Bm[b]; - for (int a = 0; a < eNoN; a++) { + // Geometric stiffness: S*grad(N_b) + const Eigen::Vector3d SNx = S * Nxm.col(b); - // Geometric stiffness - NxSNx = Nx(0,a)*S(0,0)*Nx(0,b) + Nx(1,a)*S(1,0)*Nx(0,b) + - Nx(2,a)*S(2,0)*Nx(0,b) + Nx(0,a)*S(0,1)*Nx(1,b) + - Nx(1,a)*S(1,1)*Nx(1,b) + Nx(2,a)*S(2,1)*Nx(1,b) + - Nx(0,a)*S(0,2)*Nx(2,b) + Nx(1,a)*S(1,2)*Nx(2,b) + - Nx(2,a)*S(2,2)*Nx(2,b); + for (int a = 0; a < eNoN; a++) { + NxSNx = Nxm.col(a).dot(SNx); T1 = amd*N(a)*N(b) + afu*NxSNx; // dM1/du1 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + - Bm(2,0,a)*DBm(2,0) + Bm(3,0,a)*DBm(3,0) + - Bm(4,0,a)*DBm(4,0) + Bm(5,0,a)*DBm(5,0); - - lK(0,a,b) = lK(0,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(0,a,b)) + afv*Kvis_v(0,a,b) ); + BmDBm = Bm[a].col(0).dot(DBm.col(0)); + lK(0,a,b) += w*( T1 + afu*(BmDBm + visc.du(0,a,b)) + afv*visc.dv(0,a,b) ); // dM1/du2 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,0,a)*DBm(0,1) + Bm(1,0,a)*DBm(1,1) + - Bm(2,0,a)*DBm(2,1) + Bm(3,0,a)*DBm(3,1) + - Bm(4,0,a)*DBm(4,1) + Bm(5,0,a)*DBm(5,1); - - - lK(1,a,b) = lK(1,a,b) + w*( afu*(BmDBm + Kvis_u(1,a,b)) + afv*(Kvis_v(1,a,b)) ); + BmDBm = Bm[a].col(0).dot(DBm.col(1)); + lK(1,a,b) += w*( afu*(BmDBm + visc.du(1,a,b)) + afv*visc.dv(1,a,b) ); // dM1/du3 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,0,a)*DBm(0,2) + Bm(1,0,a)*DBm(1,2) + - Bm(2,0,a)*DBm(2,2) + Bm(3,0,a)*DBm(3,2) + - Bm(4,0,a)*DBm(4,2) + Bm(5,0,a)*DBm(5,2); - - lK(2,a,b) = lK(2,a,b) + w*( afu*(BmDBm + Kvis_u(2,a,b)) + afv*Kvis_v(2,a,b) ); + BmDBm = Bm[a].col(0).dot(DBm.col(2)); + lK(2,a,b) += w*( afu*(BmDBm + visc.du(2,a,b)) + afv*visc.dv(2,a,b) ); // dM2/du1 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,1,a)*DBm(0,0) + Bm(1,1,a)*DBm(1,0) + - Bm(2,1,a)*DBm(2,0) + Bm(3,1,a)*DBm(3,0) + - Bm(4,1,a)*DBm(4,0) + Bm(5,1,a)*DBm(5,0); - - lK(dof+0,a,b) = lK(dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(3,a,b)) + afv*Kvis_v(3,a,b) ); + BmDBm = Bm[a].col(1).dot(DBm.col(0)); + lK(dof+0,a,b) += w*( afu*(BmDBm + visc.du(3,a,b)) + afv*visc.dv(3,a,b) ); // dM2/du2 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + - Bm(2,1,a)*DBm(2,1) + Bm(3,1,a)*DBm(3,1) + - Bm(4,1,a)*DBm(4,1) + Bm(5,1,a)*DBm(5,1); - - lK(dof+1,a,b) = lK(dof+1,a,b) + w*(T1 + afu*(BmDBm + Kvis_u(4,a,b)) + afv*Kvis_v(4,a,b) ); + BmDBm = Bm[a].col(1).dot(DBm.col(1)); + lK(dof+1,a,b) += w*(T1 + afu*(BmDBm + visc.du(4,a,b)) + afv*visc.dv(4,a,b) ); // dM2/du3 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,1,a)*DBm(0,2) + Bm(1,1,a)*DBm(1,2) + - Bm(2,1,a)*DBm(2,2) + Bm(3,1,a)*DBm(3,2) + - Bm(4,1,a)*DBm(4,2) + Bm(5,1,a)*DBm(5,2); - - lK(dof+2,a,b) = lK(dof+2,a,b) + w*( afu*(BmDBm + Kvis_u(5,a,b)) + afv*Kvis_v(5,a,b) ); + BmDBm = Bm[a].col(1).dot(DBm.col(2)); + lK(dof+2,a,b) += w*( afu*(BmDBm + visc.du(5,a,b)) + afv*visc.dv(5,a,b) ); // dM3/du1 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,2,a)*DBm(0,0) + Bm(1,2,a)*DBm(1,0) + - Bm(2,2,a)*DBm(2,0) + Bm(3,2,a)*DBm(3,0) + - Bm(4,2,a)*DBm(4,0) + Bm(5,2,a)*DBm(5,0); - - lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(6,a,b)) + afv*Kvis_v(6,a,b) ); + BmDBm = Bm[a].col(2).dot(DBm.col(0)); + lK(2*dof+0,a,b) += w*( afu*(BmDBm + visc.du(6,a,b)) + afv*visc.dv(6,a,b) ); // dM3/du2 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,2,a)*DBm(0,1) + Bm(1,2,a)*DBm(1,1) + - Bm(2,2,a)*DBm(2,1) + Bm(3,2,a)*DBm(3,1) + - Bm(4,2,a)*DBm(4,1) + Bm(5,2,a)*DBm(5,1); - - lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*( afu*(BmDBm + Kvis_u(7,a,b)) + afv*Kvis_v(7,a,b) ); + BmDBm = Bm[a].col(2).dot(DBm.col(1)); + lK(2*dof+1,a,b) += w*( afu*(BmDBm + visc.du(7,a,b)) + afv*visc.dv(7,a,b) ); // dM3/du3 - // Material stiffness: Bt*D*B - BmDBm = Bm(0,2,a)*DBm(0,2) + Bm(1,2,a)*DBm(1,2) + - Bm(2,2,a)*DBm(2,2) + Bm(3,2,a)*DBm(3,2) + - Bm(4,2,a)*DBm(4,2) + Bm(5,2,a)*DBm(5,2); - - lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(8,a,b)) + afv*Kvis_v(8,a,b) ); + BmDBm = Bm[a].col(2).dot(DBm.col(2)); + lK(2*dof+2,a,b) += w*( T1 + afu*(BmDBm + visc.du(8,a,b)) + afv*visc.dv(8,a,b) ); } } } diff --git a/Code/Source/solver/sv_struct.h b/Code/Source/solver/sv_struct.h index a0018305c..89d082d71 100644 --- a/Code/Source/solver/sv_struct.h +++ b/Code/Source/solver/sv_struct.h @@ -26,7 +26,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &fN, const Array &pS0l, Vector &pSl, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK); + Array &lR, Array3 &lK, const bool recompute_visc); void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const double w, const Vector &N, const Array &Nx, @@ -35,7 +35,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &fN, const Array &pS0l, Vector &pSl, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK); + Array &lR, Array3 &lK, const bool recompute_visc); }; #endif diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 87f5efa4d..a7f3c0908 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -27,7 +27,7 @@ #include "mat_models.h" #include "nn.h" #include "utils.h" - +#include #include namespace ustruct { @@ -268,6 +268,14 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const ya_l_s = 0.0; ya_l_n = 0.0; + if (lM.fN.size() != 0) { + for (int iFn = 0; iFn < nFn; iFn++) { + for (int i = 0; i < nsd; i++) { + fN(i,iFn) = lM.fN(i+nsd*iFn,e); + } + } + } + for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); ptr(a) = Ac; @@ -283,14 +291,6 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const yl(i,a) = Yg(i,Ac); } - if (lM.fN.size() != 0) { - for (int iFn = 0; iFn < nFn; iFn++) { - for (int i = 0; i < nsd; i++) { - fN(i,iFn) = lM.fN(i+nsd*iFn,e); - } - } - } - if (eq.dmn[cDmn].active_stress != nullptr) { ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; @@ -327,7 +327,11 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array ksix(nsd,nsd); for (int g = 0; g < fs[0].nG; g++) { - if (g == 0 || !fs[0].lShpF) { + // Shape function gradients and the viscous response are constant + // within linear triangles and tetrahedra. triangles). + const bool recompute_visc = (g == 0 || !fs[0].lShpF); + + if (recompute_visc) { auto Nx = fs[0].Nx.slice(g); nn::gnn(fs[0].eNoN, nsd, nsd, Nx, xwl, Nwx, Jac, ksix); if (utils::is_zero(Jac)) { @@ -342,14 +346,14 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const auto N1 = fs[1].N.col(g); ustruct_3d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd); + ya_l_n, lR, lK, lKd, recompute_visc); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_2d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd); + ya_l_n, lR, lK, lKd, recompute_visc); } } // for g = 0 to fs[0].nG @@ -448,19 +452,16 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - Vector fb(2); - fb[0] = dmn.prop[PhysicalPropertyType::f_x]; - fb[1] = dmn.prop[PhysicalPropertyType::f_y]; - fb[2] = dmn.prop[PhysicalPropertyType::f_z]; + const Eigen::Vector2d fb{dmn.prop[PhysicalPropertyType::f_x], + dmn.prop[PhysicalPropertyType::f_y]}; double am = eq.am; double af = eq.af * eq.gam * dt; double afm = af / am; - // {i,j} := velocity dofs; {k} := pressure dof + // Velocity dofs start at i; k is the pressure dof. int i = eq.s; - int j = i + 1; - int k = j + 1; + int k = i + 2; #ifdef debug_ustruct_2d_c dmsg << "am: " << am; @@ -469,35 +470,24 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in dmsg << "i: " << i; #endif - // Inertia (velocity and acceleration), body force, fiber directions, - // and deformation tensor (F) at integration point - // - Vector vd{-fb[0], -fb[1]}; - Vector v(2); - Array vx(2,2), F(2,2); - F(0,0) = 1.0; - F(1,1) = 1.0; + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nwxm = eigen_view<2>(Nwx); // grad(N_a) per column + const auto Nwm = eigen_view(Nw); // shape functions + const auto disp = eigen_view(dl).middleRows<2>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<2>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<2>(i); // nodal accelerations + const auto bfm = eigen_view<2>(bfl); // nodal body force - for (int a = 0; a < eNoNw; a++) { - v(0) = v(0) + Nw(a)*yl(i,a); - v(1) = v(1) + Nw(a)*yl(j,a); + // Velocity and inertia at this Gauss point + const Eigen::Vector2d v = vel * Nwm; + const Eigen::Vector2d vd = (acc - bfm) * Nwm - fb; - vd(0) = vd(0) + Nw(a)*(al(i,a)-bfl(0,a)); - vd(1) = vd(1) + Nw(a)*(al(j,a)-bfl(1,a)); + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<2> vx = vel * Nwxm.transpose(); + const Matrix<2> F = Matrix<2>::Identity() + disp * Nwxm.transpose(); - vx(0,0) = vx(0,0) + Nwx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nwx(1,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nwx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nwx(1,a)*yl(j,a); - - F(0,0) = F(0,0) + Nwx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nwx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nwx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nwx(1,a)*dl(j,a); - } - - double Jac = mat_fun::mat_det(F, 2); - auto Fi = mat_fun::mat_inv(F, 2); + double Jac = F.determinant(); + const Matrix<2> Fi = F.inverse(); // Pressure and its gradients // @@ -508,8 +498,8 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int a = 0; a < eNoNq; a++) { p = p + Nq(a)*yl(k,a); pd = pd + Nq(a)*al(k,a); - px(0) = px(0) + Nqx(0,a)*yl(k,a); - px(1) = px(1) + Nqx(1,a)*yl(k,a); + px(0) += Nqx(0,a)*yl(k,a); + px(1) += Nqx(1,a)*yl(k,a); } // Compute rho and beta depending on the volumetric penalty model @@ -532,30 +522,13 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in tauC = 0.0; } - Array NwxFi(2,eNoNw); - - for (int a = 0; a < eNoNw; a++) { - NwxFi(0,a) = Nwx(0,a)*Fi(0,0) + Nwx(1,a)*Fi(1,0); - NwxFi(1,a) = Nwx(0,a)*Fi(0,1) + Nwx(1,a)*Fi(1,1); - } - - Array NqxFi(2,eNoNw); - - for (int a = 0; a < eNoNq; a++) { - NqxFi(0,a) = Nqx(0,a)*Fi(0,0) + Nqx(1,a)*Fi(1,0); - NqxFi(1,a) = Nqx(0,a)*Fi(0,1) + Nqx(1,a)*Fi(1,1); - } + const NodalMatrix<2> NwxFi = Fi.transpose() * Nwxm; - Array VxFi(2,2); + const NodalMatrix<2> NqxFi = Fi.transpose() * eigen_view<2>(Nqx); - VxFi(0,0) = vx(0,0)*Fi(0,0) + vx(0,1)*Fi(1,0); - VxFi(0,1) = vx(0,0)*Fi(0,1) + vx(0,1)*Fi(1,1); - VxFi(1,0) = vx(1,0)*Fi(0,0) + vx(1,1)*Fi(1,0); - VxFi(1,1) = vx(1,0)*Fi(0,1) + vx(1,1)*Fi(1,1); + const Matrix<2> VxFi = vx * Fi; - Vector PxFi(2); - PxFi(0) = px(0)*Fi(0,0) + px(1)*Fi(1,0); - PxFi(1) = px(0)*Fi(0,1) + px(1)*Fi(1,1); + const Eigen::Vector2d PxFi = Fi.transpose() * eigen_view<2>(px); double rC = beta*pd + VxFi(0,0) + VxFi(1,1); @@ -565,21 +538,15 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // Local residual // - Vector rMNqx(eNoNq); + const NodalVector rMNqx = NqxFi.transpose() * eigen_view<2>(rM); for (int a = 0; a < eNoNq; a++) { - rMNqx(a) = rM(0)*NqxFi(0,a) + rM(1)*NqxFi(1,a); - lR(2,a) = lR(2,a) + w*Jac*(Nq(a)*rC + tauM*rMNqx(a)); + lR(2,a) += w*Jac*(Nq(a)*rC + tauM*rMNqx(a)); } - Vector rMNwx(eNoNw); - Array VxNwx(3,eNoNw); + const NodalMatrix<2> VxNwx = VxFi.transpose() * NwxFi; - for (int a = 0; a < eNoNw; a++) { - rMNwx(a) = rM(0)*NwxFi(0,a) + rM(1)*NwxFi(1,a); - VxNwx(0,a) = VxFi(0,0)*NwxFi(0,a) + VxFi(1,0)*NwxFi(1,a); - VxNwx(1,a) = VxFi(0,1)*NwxFi(0,a) + VxFi(1,1)*NwxFi(1,a); - } + const NodalVector rMNwx = NwxFi.transpose() * eigen_view<2>(rM); // Tangent (stiffness) matrices // @@ -587,7 +554,7 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int b = 0; b < eNoNw; b++) { for (int a = 0; a < eNoNq; a++) { - NxNx = NqxFi(0,a)*NwxFi(0,b) + NqxFi(1,a)*NwxFi(1,b); + NxNx = NqxFi.col(a).dot(NwxFi.col(b)); // dC/dV_1 + af/am *dC/dU_1 // @@ -595,31 +562,31 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = tauM*(rMNqx(a)*NwxFi(0,b) - rMNwx(b)*NqxFi(0,a)); T2 = -tauM*NxNx*PxFi(0); Ku = w*af*Jac*(T0 + T1 + T2); - lKd(4,a,b) = lKd(4,a,b) + Ku; + lKd(4,a,b) += Ku; T1 = (am*tauM*rho)*NqxFi(0,a)*Nw(b) + af*Nq(a)*NwxFi(0,b); - lK(6,a,b) = lK(6,a,b) + w*Jac*T1 + afm*Ku; + lK(6,a,b) += w*Jac*T1 + afm*Ku; // dC/dV_2 + af/am *dC/dU_2 T0 = Nq(a)*(rC*NwxFi(1,b) - VxNwx(1,b)); T1 = tauM*(rMNqx(a)*NwxFi(1,b) - rMNwx(b)*NqxFi(1,a)); T2 = -tauM*NxNx*PxFi(1); Ku = w*af*Jac*(T0 + T1 + T2); - lKd(5,a,b) = lKd(5,a,b) + Ku; + lKd(5,a,b) += Ku; T1 = (am*tauM*rho)*NqxFi(1,a)*Nw(b) + af*Nq(a)*NwxFi(1,b); - lK(8,a,b) = lK(8,a,b) + w*Jac*T1 + afm*Ku; + lK(8,a,b) += w*Jac*T1 + afm*Ku; } } for (int b = 0; b < eNoNq; b++) { for (int a = 0; a < eNoNq; a++) { // dC/dP - NxNx = NqxFi(0,a)*NqxFi(0,b) + NqxFi(1,a)*NqxFi(1,b); + NxNx = NqxFi.col(a).dot(NqxFi.col(b)); T0 = (am*beta + af*dbeta*pd)*Nq(a)*Nq(b); - T1 = NqxFi(0,a)*vd(0) + NqxFi(1,a)*vd(1); + T1 = NqxFi.col(a).dot(vd); T2 = T0 + af*tauM*(NxNx + drho*T1*Nq(b)); - lK(9,a,b) = lK(9,a,b) + w*Jac*T2; + lK(9,a,b) += w*Jac*T2; } } } @@ -652,20 +619,17 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - Vector fb(3); - fb[0] = dmn.prop[PhysicalPropertyType::f_x]; - fb[1] = dmn.prop[PhysicalPropertyType::f_y]; - fb[2] = dmn.prop[PhysicalPropertyType::f_z]; + const Eigen::Vector3d fb{dmn.prop[PhysicalPropertyType::f_x], + dmn.prop[PhysicalPropertyType::f_y], + dmn.prop[PhysicalPropertyType::f_z]}; double am = eq.am; double af = eq.af * eq.gam * dt; double afm = af / am; - // {i,j} := velocity dofs; {k} := pressure dof + // Velocity dofs start at i; l is the pressure dof. int i = eq.s; - int j = i + 1; - int k = j + 1; - int l = k + 1; + int l = i + 3; #ifdef debug_ustruct_3d_c dmsg << "am: " << am; @@ -674,52 +638,24 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in dmsg << "i: " << i; #endif - // Inertia (velocity and acceleration), body force, fiber directions, - // and deformation tensor (F) at integration point - // - Vector vd{-fb[0], -fb[1], -fb[2]}; - Vector v(3); - Array vx(3,3), F(3,3); - F(0,0) = 1.0; - F(1,1) = 1.0; - F(2,2) = 1.0; - - for (int a = 0; a < eNoNw; a++) { - v(0) = v(0) + Nw(a)*yl(i,a); - v(1) = v(1) + Nw(a)*yl(j,a); - v(2) = v(2) + Nw(a)*yl(k,a); + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nwxm = eigen_view<3>(Nwx); // grad(N_a) per column + const auto Nwm = eigen_view(Nw); // shape functions + const auto disp = eigen_view(dl).middleRows<3>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<3>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<3>(i); // nodal accelerations + const auto bfm = eigen_view<3>(bfl); // nodal body force - vd(0) = vd(0) + Nw(a)*(al(i,a)-bfl(0,a)); - vd(1) = vd(1) + Nw(a)*(al(j,a)-bfl(1,a)); - vd(2) = vd(2) + Nw(a)*(al(k,a)-bfl(2,a)); + // Velocity and inertia at this Gauss point + const Eigen::Vector3d v = vel * Nwm; + const Eigen::Vector3d vd = (acc - bfm) * Nwm - fb; - vx(0,0) = vx(0,0) + Nwx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nwx(1,a)*yl(i,a); - vx(0,2) = vx(0,2) + Nwx(2,a)*yl(i,a); + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<3> vx = vel * Nwxm.transpose(); + const Matrix<3> F = Matrix<3>::Identity() + disp * Nwxm.transpose(); - vx(1,0) = vx(1,0) + Nwx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nwx(1,a)*yl(j,a); - vx(1,2) = vx(1,2) + Nwx(2,a)*yl(j,a); - - vx(2,0) = vx(2,0) + Nwx(0,a)*yl(k,a); - vx(2,1) = vx(2,1) + Nwx(1,a)*yl(k,a); - vx(2,2) = vx(2,2) + Nwx(2,a)*yl(k,a); - - F(0,0) = F(0,0) + Nwx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nwx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nwx(2,a)*dl(i,a); - - F(1,0) = F(1,0) + Nwx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nwx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nwx(2,a)*dl(j,a); - - F(2,0) = F(2,0) + Nwx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nwx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nwx(2,a)*dl(k,a); - } - - double Jac = mat_fun::mat_det(F, 3); - auto Fi = mat_fun::mat_inv(F, 3); + double Jac = F.determinant(); + const Matrix<3> Fi = F.inverse(); // Pressure and its gradients // @@ -730,9 +666,9 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int a = 0; a < eNoNq; a++) { p = p + Nq(a)*yl(l,a); pd = pd + Nq(a)*al(l,a); - px(0) = px(0) + Nqx(0,a)*yl(l,a); - px(1) = px(1) + Nqx(1,a)*yl(l,a); - px(2) = px(2) + Nqx(2,a)*yl(l,a); + px(0) += Nqx(0,a)*yl(l,a); + px(1) += Nqx(1,a)*yl(l,a); + px(2) += Nqx(2,a)*yl(l,a); } // Compute rho and beta depending on the volumetric penalty model @@ -755,40 +691,13 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in tauC = 0.0; } - Array NwxFi(3,eNoNw); + const NodalMatrix<3> NwxFi = Fi.transpose() * Nwxm; - for (int a = 0; a < eNoNw; a++) { - NwxFi(0,a) = Nwx(0,a)*Fi(0,0) + Nwx(1,a)*Fi(1,0) + Nwx(2,a)*Fi(2,0); - NwxFi(1,a) = Nwx(0,a)*Fi(0,1) + Nwx(1,a)*Fi(1,1) + Nwx(2,a)*Fi(2,1); - NwxFi(2,a) = Nwx(0,a)*Fi(0,2) + Nwx(1,a)*Fi(1,2) + Nwx(2,a)*Fi(2,2); - } + const NodalMatrix<3> NqxFi = Fi.transpose() * eigen_view<3>(Nqx); - Array NqxFi(3,eNoNw); + const Matrix<3> VxFi = vx * Fi; - for (int a = 0; a < eNoNq; a++) { - NqxFi(0,a) = Nqx(0,a)*Fi(0,0) + Nqx(1,a)*Fi(1,0) + Nqx(2,a)*Fi(2,0); - NqxFi(1,a) = Nqx(0,a)*Fi(0,1) + Nqx(1,a)*Fi(1,1) + Nqx(2,a)*Fi(2,1); - NqxFi(2,a) = Nqx(0,a)*Fi(0,2) + Nqx(1,a)*Fi(1,2) + Nqx(2,a)*Fi(2,2); - } - - Array VxFi(3,3); - - VxFi(0,0) = vx(0,0)*Fi(0,0) + vx(0,1)*Fi(1,0) + vx(0,2)*Fi(2,0); - VxFi(0,1) = vx(0,0)*Fi(0,1) + vx(0,1)*Fi(1,1) + vx(0,2)*Fi(2,1); - VxFi(0,2) = vx(0,0)*Fi(0,2) + vx(0,1)*Fi(1,2) + vx(0,2)*Fi(2,2); - - VxFi(1,0) = vx(1,0)*Fi(0,0) + vx(1,1)*Fi(1,0) + vx(1,2)*Fi(2,0); - VxFi(1,1) = vx(1,0)*Fi(0,1) + vx(1,1)*Fi(1,1) + vx(1,2)*Fi(2,1); - VxFi(1,2) = vx(1,0)*Fi(0,2) + vx(1,1)*Fi(1,2) + vx(1,2)*Fi(2,2); - - VxFi(2,0) = vx(2,0)*Fi(0,0) + vx(2,1)*Fi(1,0) + vx(2,2)*Fi(2,0); - VxFi(2,1) = vx(2,0)*Fi(0,1) + vx(2,1)*Fi(1,1) + vx(2,2)*Fi(2,1); - VxFi(2,2) = vx(2,0)*Fi(0,2) + vx(2,1)*Fi(1,2) + vx(2,2)*Fi(2,2); - - Vector PxFi(3); - PxFi(0) = px(0)*Fi(0,0) + px(1)*Fi(1,0) + px(2)*Fi(2,0); - PxFi(1) = px(0)*Fi(0,1) + px(1)*Fi(1,1) + px(2)*Fi(2,1); - PxFi(2) = px(0)*Fi(0,2) + px(1)*Fi(1,2) + px(2)*Fi(2,2); + const Eigen::Vector3d PxFi = Fi.transpose() * eigen_view<3>(px); double rC = beta*pd + VxFi(0,0) + VxFi(1,1) + VxFi(2,2); @@ -799,22 +708,15 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // Local residual // - Vector rMNqx(eNoNq); + const NodalVector rMNqx = NqxFi.transpose() * eigen_view<3>(rM); for (int a = 0; a < eNoNq; a++) { - rMNqx(a) = rM(0)*NqxFi(0,a) + rM(1)*NqxFi(1,a) + rM(2)*NqxFi(2,a); - lR(3,a) = lR(3,a) + w*Jac*(Nq(a)*rC + tauM*rMNqx(a)); + lR(3,a) += w*Jac*(Nq(a)*rC + tauM*rMNqx(a)); } - Vector rMNwx(eNoNw); - Array VxNwx(3,eNoNw); + const NodalMatrix<3> VxNwx = VxFi.transpose() * NwxFi; - for (int a = 0; a < eNoNw; a++) { - rMNwx(a) = rM(0)*NwxFi(0,a) + rM(1)*NwxFi(1,a) + rM(2)*NwxFi(2,a); - VxNwx(0,a) = VxFi(0,0)*NwxFi(0,a) + VxFi(1,0)*NwxFi(1,a) + VxFi(2,0)*NwxFi(2,a); - VxNwx(1,a) = VxFi(0,1)*NwxFi(0,a) + VxFi(1,1)*NwxFi(1,a) + VxFi(2,1)*NwxFi(2,a); - VxNwx(2,a) = VxFi(0,2)*NwxFi(0,a) + VxFi(1,2)*NwxFi(1,a) + VxFi(2,2)*NwxFi(2,a); - } + const NodalVector rMNwx = NwxFi.transpose() * eigen_view<3>(rM); // Tangent (stiffness) matrices // @@ -822,7 +724,7 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int b = 0; b < eNoNw; b++) { for (int a = 0; a < eNoNq; a++) { - NxNx = NqxFi(0,a)*NwxFi(0,b) + NqxFi(1,a)*NwxFi(1,b) + NqxFi(2,a)*NwxFi(2,b); + NxNx = NqxFi.col(a).dot(NwxFi.col(b)); // dC/dV_1 + af/am *dC/dU_1 // @@ -830,20 +732,20 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = tauM*(rMNqx(a)*NwxFi(0,b) - rMNwx(b)*NqxFi(0,a)); T2 = -tauM*NxNx*PxFi(0); Ku = w*af*Jac*(T0 + T1 + T2); - lKd(9,a,b) = lKd(9,a,b) + Ku; + lKd(9,a,b) += Ku; T1 = (am*tauM*rho)*NqxFi(0,a)*Nw(b) + af*Nq(a)*NwxFi(0,b); - lK(12,a,b) = lK(12,a,b) + w*Jac*T1 + afm*Ku; + lK(12,a,b) += w*Jac*T1 + afm*Ku; // dC/dV_2 + af/am *dC/dU_2 T0 = Nq(a)*(rC*NwxFi(1,b) - VxNwx(1,b)); T1 = tauM*(rMNqx(a)*NwxFi(1,b) - rMNwx(b)*NqxFi(1,a)); T2 = -tauM*NxNx*PxFi(1); Ku = w*af*Jac*(T0 + T1 + T2); - lKd(10,a,b) = lKd(10,a,b) + Ku; + lKd(10,a,b) += Ku; T1 = (am*tauM*rho)*NqxFi(1,a)*Nw(b) + af*Nq(a)*NwxFi(1,b); - lK(13,a,b) = lK(13,a,b) + w*Jac*T1 + afm*Ku; + lK(13,a,b) += w*Jac*T1 + afm*Ku; // dC/dV_3 + af/am *dC/dU_3 // @@ -851,21 +753,21 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = tauM*(rMNqx(a)*NwxFi(2,b) - rMNwx(b)*NqxFi(2,a)); T2 = -tauM*NxNx*PxFi(2); Ku = w*af*Jac*(T0 + T1 + T2); - lKd(11,a,b) = lKd(11,a,b) + Ku; + lKd(11,a,b) += Ku; T1 = (am*tauM*rho)*NqxFi(2,a)*Nw(b) + af*Nq(a)*NwxFi(2,b); - lK(14,a,b) = lK(14,a,b) + w*Jac*T1 + afm*Ku; + lK(14,a,b) += w*Jac*T1 + afm*Ku; } } for (int b = 0; b < eNoNq; b++) { for (int a = 0; a < eNoNq; a++) { // dC/dP - NxNx = NqxFi(0,a)*NqxFi(0,b) + NqxFi(1,a)*NqxFi(1,b) + NqxFi(2,a)*NqxFi(2,b); + NxNx = NqxFi.col(a).dot(NqxFi.col(b)); T0 = (am*beta + af*dbeta*pd)*Nq(a)*Nq(b); - T1 = NqxFi(0,a)*vd(0) + NqxFi(1,a)*vd(1) + NqxFi(2,a)*vd(2); + T1 = NqxFi.col(a).dot(vd); T2 = T0 + af*tauM*(NxNx + drho*T1*Nq(b)); - lK(15,a,b) = lK(15,a,b) + w*Jac*T2; + lK(15,a,b) += w*Jac*T2; } } } @@ -880,7 +782,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &dl, const Array &bfl, const Array &fN, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK, Array3 &lKd) { + Array &lR, Array3 &lK, Array3 &lKd, + const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -901,18 +804,16 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Define parameters // - Vector fb(2); - fb[0] = dmn.prop[PhysicalPropertyType::f_x]; - fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + const Eigen::Vector2d fb{dmn.prop[PhysicalPropertyType::f_x], + dmn.prop[PhysicalPropertyType::f_y]}; double am = eq.am; double af = eq.af * eq.gam * dt; double afm = af / am; - // {i,j} := velocity dofs; {k} := pressure dof + // Velocity dofs start at i; k is the pressure dof. int i = eq.s; - int j = i + 1; - int k = j + 1; + int k = i + 2; #ifdef debug_ustruct_2d_m dmsg << "am: " << am; @@ -921,44 +822,30 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, dmsg << "i: " << i; #endif - // Inertia (velocity and acceleration), body force, fiber directions, - // and deformation tensor (F) at integration point - // - Vector vd{-fb[0], -fb[1]}; - Vector v(2); - Array vx(2,2), F(2,2); - - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - - F(0,0) = 1.0; - F(1,1) = 1.0; - - for (int a = 0; a < eNoNw; a++) { - v(0) = v(0) + Nw(a)*yl(i,a); - v(1) = v(1) + Nw(a)*yl(j,a); - - vd(0) = vd(0) + Nw(a)*(al(i,a)-bfl(0,a)); - vd(1) = vd(1) + Nw(a)*(al(j,a)-bfl(1,a)); + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nwxm = eigen_view<2>(Nwx); // grad(N_a) per column + const auto Nwm = eigen_view(Nw); // shape functions + const auto disp = eigen_view(dl).middleRows<2>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<2>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<2>(i); // nodal accelerations + const auto bfm = eigen_view<2>(bfl); // nodal body force + auto lRv = eigen_view_mutable(lR).topRows<2>(); // rows this kernel adds to - vx(0,0) = vx(0,0) + Nwx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nwx(1,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nwx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nwx(1,a)*yl(j,a); + // Velocity and inertia at this Gauss point + const Eigen::Vector2d v = vel * Nwm; + const Eigen::Vector2d vd = (acc - bfm) * Nwm - fb; - F(0,0) = F(0,0) + Nwx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nwx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nwx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nwx(1,a)*dl(j,a); + // Active stress activation along fiber, sheet and sheet-normal + const double ya_g_f = eigen_view(ya_l_f).dot(Nwm); + const double ya_g_s = eigen_view(ya_l_s).dot(Nwm); + const double ya_g_n = eigen_view(ya_l_n).dot(Nwm); - ya_g_f = ya_g_f + Nw(a) * ya_l_f(a); - ya_g_s = ya_g_s + Nw(a) * ya_l_s(a); - ya_g_n = ya_g_n + Nw(a) * ya_l_n(a); - } + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<2> vx = vel * Nwxm.transpose(); + const Matrix<2> F = Matrix<2>::Identity() + disp * Nwxm.transpose(); - double Jac = mat_fun::mat_det(F, 2); - auto Fi = mat_fun::mat_inv(F, 2); + double Jac = F.determinant(); + const Matrix<2> Fi = F.inverse(); // Pressure and its time derivative // @@ -972,17 +859,18 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Compute deviatoric 2nd Piola-Kirchhoff stress tensor (Siso) and // isochoric elasticity tensor in Voigt notation (Dm) - Array Siso(2,2), Dm(3,3); + Matrix<2> Siso; + Matrix<3> Dm; double Ja = 0; - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g_f, + mat_models::compute_pk2cc<2>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<2>(fN), ya_g_f, ya_g_s, ya_g_n, Siso, Dm, Ja); - // Viscous 2nd Piola-Kirchhoff stress and tangent contributions - Array Svis(2,2); - Array3 Kvis_u(4, eNoNw, eNoNw); - Array3 Kvis_v(4, eNoNw, eNoNw); - - mat_models::compute_visc_stress_and_tangent(dmn, eNoNw, Nwx, vx, F, Svis, Kvis_u, Kvis_v); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. Reuse the + // previous Gauss point's when shape function gradients are constant within an + // element (e.g. linear triangles and tetrahedra). + static mat_models::ViscousResponse<2> visc; + + visc.update(dmn, eNoNw, Nwx, vx, F, recompute_visc); // Compute rho and beta depending on the volumetric penalty model // @@ -1004,62 +892,41 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, } // Total isochoric 2nd Piola-Kirchhoff stress (Elastic + Viscous) - Siso = Siso + Svis; + Siso += visc.S(); // Deviatoric 1st Piola-Kirchhoff tensor (P) // - auto Pdev = mat_fun::mat_mul(F, Siso); - + const Matrix<2> Pdev = F * Siso; // Shape function gradients in the current configuration // - Array NxFi(2,eNoNw); - - for (int a = 0; a < eNoNw; a++) { - NxFi(0,a) = Nwx(0,a)*Fi(0,0) + Nwx(1,a)*Fi(1,0); - NxFi(1,a) = Nwx(0,a)*Fi(0,1) + Nwx(1,a)*Fi(1,1); - } + const NodalMatrix<2> NxFi = Fi.transpose() * Nwxm; // Velocity gradient in current configuration - auto VxFi = mat_mul(vx, Fi); - double rC = beta*pd + VxFi(1,1) + VxFi(2,2); + const Matrix<2> VxFi = vx * Fi; + double rC = beta*pd + VxFi(0,0) + VxFi(1,1); double rCl = -p + tauC*rC; // Local residual // - for (int a = 0; a < eNoNw; a++) { - double T1 = Jac*rho*vd(0)*Nw(a); - double T2 = Pdev(0,0)*Nwx(0,a) + Pdev(0,1)*Nwx(1,a); - double T3 = Jac*rCl*NxFi(0,a); - lR(0,a) = lR(0,a) + w*(T1 + T2 + T3); - - T1 = Jac*rho*vd(1)*Nw(a); - T2 = Pdev(1,0)*Nwx(0,a) + Pdev(1,1)*Nwx(1,a); - T3 = Jac*rCl*NxFi(1,a); - lR(1,a) = lR(1,a) + w*(T1 + T2 + T3); - } + // Inertia, the divergence of Pdev, and the pressure/volumetric term + lRv += w * (Jac*rho * vd * Nwm.transpose() + Pdev * Nwxm + Jac*rCl * NxFi); - // Auxilary quantities for computing stiffness tensors + // Strain-displacement matrix; Bm[a] maps node a to Voigt strain // - Array3 Bm(3,2,eNoNw); + std::array, consts::maxNoN> Bm; - for (int a = 0; a < eNoNw; a++) { - Bm(0,0,a) = Nwx(0,a)*F(0,0); - Bm(0,1,a) = Nwx(0,a)*F(1,0); + const Matrix<2> Ft = F.transpose(); - Bm(1,0,a) = Nwx(1,a)*F(0,1); - Bm(1,1,a) = Nwx(1,a)*F(1,1); + for (int a = 0; a < eNoNw; a++) { + const auto g = Nwxm.col(a); // grad(N_a) - Bm(2,0,a) = Nwx(2,a)*F(0,2) + F(0,0)*Nwx(1,a); - Bm(2,1,a) = Nwx(2,a)*F(1,2) + F(1,0)*Nwx(1,a); + Bm[a].row(0) = g(0) * Ft.row(0); // dE_11 + Bm[a].row(1) = g(1) * Ft.row(1); // dE_22 + Bm[a].row(2) = g(0) * Ft.row(1) + g(1) * Ft.row(0); // 2 dE_12 } - Array VxNx(2,eNoNw); - - for (int a = 0; a < eNoNw; a++) { - VxNx(0,a) = VxFi(0,0)*NxFi(0,a) + VxFi(1,0)*NxFi(1,a); - VxNx(1,a) = VxFi(0,1)*NxFi(0,a) + VxFi(1,1)*NxFi(1,a); - } + const NodalMatrix<2> VxNx = VxFi.transpose() * NxFi; // Tangent (stiffness) matrices // @@ -1067,80 +934,71 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1{0.0}, T2{0.0}, T3{0.0}, Tv{0.0}, Ku{0.0}; - Array DBm(3,2); - for (int b = 0; b < eNoNw; b++) { - for (int a = 0; a < eNoNw; a++) { - NxSNx = Nwx(0,a)*Siso(0,0)*Nwx(0,b) - + Nwx(0,a)*Siso(0,1)*Nwx(1,b) - + Nwx(1,a)*Siso(1,0)*Nwx(0,b) - + Nwx(1,a)*Siso(1,1)*Nwx(1,b); - - DBm(0,0) = Dm(0,0)*Bm(0,0,b) + Dm(0,1)*Bm(1,0,b) + Dm(0,2)*Bm(2,0,b); - DBm(0,1) = Dm(0,0)*Bm(0,1,b) + Dm(0,1)*Bm(1,1,b) + Dm(0,2)*Bm(2,1,b); - DBm(1,0) = Dm(1,0)*Bm(0,0,b) + Dm(1,1)*Bm(1,0,b) + Dm(1,2)*Bm(2,0,b); - DBm(1,1) = Dm(1,0)*Bm(0,1,b) + Dm(1,1)*Bm(1,1,b) + Dm(1,2)*Bm(2,1,b); + const Eigen::Matrix DBm = Dm * Bm[b]; - DBm(2,0) = Dm(2,0)*Bm(0,0,b) + Dm(2,1)*Bm(1,0,b) + Dm(2,2)*Bm(2,0,b); - DBm(2,1) = Dm(2,0)*Bm(0,1,b) + Dm(2,1)*Bm(1,1,b) + Dm(2,2)*Bm(2,1,b); + // Geometric stiffness: Siso*grad(N_b) + const Eigen::Vector2d SisoNx = Siso * Nwxm.col(b); + for (int a = 0; a < eNoNw; a++) { + NxSNx = Nwxm.col(a).dot(SisoNx); // dM1_dV1 + af/am *dM_1/dU_1 // - BtDB = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + Bm(2,0,a)*DBm(2,0); + BtDB = Bm[a].col(0).dot(DBm.col(0)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(0,b); - Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(0,a,b)); - lKd(0,a,b) = lKd(0,a,b) + Ku; + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(0,a,b)); + lKd(0,a,b) += Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(0,a)*NxFi(0,b); - Tv = af*Kvis_v(0,a,b); - lK(0,a,b) = lK(0,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(0,a,b); + lK(0,a,b) += w*(T2 + Tv) + afm*Ku; // dM_1/dV_2 + af/am *dM_1/dU_2 // - BtDB = Bm(0,0,a)*DBm(0,1) + Bm(1,0,a)*DBm(1,1) + Bm(2,0,a)*DBm(2,1); + BtDB = Bm[a].col(0).dot(DBm.col(1)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(1,b) - NxFi(1,a)*NxFi(0,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(1,a,b)); - lKd(1,a,b) = lKd(1,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(1,a,b)); + lKd(1,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(1,b); - Tv = af*Kvis_v(1,a,b); - lK(1,a,b) = lK(1,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(1,a,b); + lK(1,a,b) += w*(T2 + Tv) + afm*Ku; // dM_2/dV_1 + af/am *dM_2/dU_1 // - BtDB = Bm(0,1,a)*DBm(0,0) + Bm(1,1,a)*DBm(1,0) + Bm(2,1,a)*DBm(2,0); + BtDB = Bm[a].col(1).dot(DBm.col(0)); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(0,b) - NxFi(0,a)*NxFi(1,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(2,a,b)); - lKd(2,a,b) = lKd(2,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(2,a,b)); + lKd(2,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(0,b); - Tv = af*Kvis_v(2,a,b); - lK(3,a,b) = lK(3,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(2,a,b); + lK(3,a,b) += w*(T2 + Tv) + afm*Ku; // dM_2/dV_2 + af/am *dM_2/dU_2 // - BtDB = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1); + BtDB = Bm[a].col(1).dot(DBm.col(1)); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(1,b); - Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(3,a,b)); - lKd(3,a,b) = lKd(3,a,b) + Ku; + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(3,a,b)); + lKd(3,a,b) += Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(1,a)*NxFi(1,b); - Tv = af*Kvis_v(3,a,b); - lK(4,a,b) = lK(4,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(3,a,b); + lK(4,a,b) += w*(T2 + Tv) + afm*Ku; } } @@ -1151,11 +1009,11 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // dM_0/dP T0 = am*tauC*beta + af*(tauC*dbeta*pd - 1.0); T1 = T0*NxFi(0,a)*Nq(b) + af*drho*vd(0)*Nw(a)*Nq(b); - lK(2,a,b) = lK(2,a,b) + w*Jac*T1; + lK(2,a,b) += w*Jac*T1; // dM_1/dP T1 = T0*NxFi(1,a)*Nq(b) + af*drho*vd(1)*Nw(a)*Nq(b); - lK(6,a,b) = lK(6,a,b) + w*Jac*T1; + lK(6,a,b) += w*Jac*T1; } } } @@ -1170,7 +1028,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &dl, const Array &bfl, const Array &fN, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK, Array3 &lKd) { + Array &lR, Array3 &lK, Array3 &lKd, + const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -1191,20 +1050,17 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Define parameters - Vector fb(3); - fb[0] = dmn.prop[PhysicalPropertyType::f_x]; - fb[1] = dmn.prop[PhysicalPropertyType::f_y]; - fb[2] = dmn.prop[PhysicalPropertyType::f_z]; + const Eigen::Vector3d fb{dmn.prop[PhysicalPropertyType::f_x], + dmn.prop[PhysicalPropertyType::f_y], + dmn.prop[PhysicalPropertyType::f_z]}; double am = eq.am; double af = eq.af * eq.gam * dt; double afm = af / am; - // {i,j} := velocity dofs; {k} := pressure dof + // Velocity dofs start at i; l is the pressure dof. int i = eq.s; - int j = i + 1; - int k = j + 1; - int l = k + 1; + int l = i + 3; #ifdef debug_ustruct_3d_m dmsg << "fb: " << fb; @@ -1214,61 +1070,30 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, dmsg << "i: " << i; #endif - // Inertia (velocity and acceleration), body force, fiber directions, - // and deformation tensor (F) at integration point - // - Vector vd{-fb[0], -fb[1], -fb[2]}; - Vector v(3); - Array vx(3,3), F(3,3); - - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - - F(0,0) = 1.0; - F(1,1) = 1.0; - F(2,2) = 1.0; - - for (int a = 0; a < eNoNw; a++) { - v(0) = v(0) + Nw(a)*yl(i,a); - v(1) = v(1) + Nw(a)*yl(j,a); - v(2) = v(2) + Nw(a)*yl(k,a); - - vd(0) = vd(0) + Nw(a)*(al(i,a)-bfl(0,a)); - vd(1) = vd(1) + Nw(a)*(al(j,a)-bfl(1,a)); - vd(2) = vd(2) + Nw(a)*(al(k,a)-bfl(2,a)); + // This element's nodal fields, as Eigen views over the caller's storage + const auto Nwxm = eigen_view<3>(Nwx); // grad(N_a) per column + const auto Nwm = eigen_view(Nw); // shape functions + const auto disp = eigen_view(dl).middleRows<3>(i); // nodal displacements + const auto vel = eigen_view(yl).middleRows<3>(i); // nodal velocities + const auto acc = eigen_view(al).middleRows<3>(i); // nodal accelerations + const auto bfm = eigen_view<3>(bfl); // nodal body force + auto lRv = eigen_view_mutable(lR).topRows<3>(); // rows this kernel adds to - vx(0,0) = vx(0,0) + Nwx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nwx(1,a)*yl(i,a); - vx(0,2) = vx(0,2) + Nwx(2,a)*yl(i,a); + // Velocity, and the inertia less body force, at this Gauss point + const Eigen::Vector3d v = vel * Nwm; + const Eigen::Vector3d vd = (acc - bfm) * Nwm - fb; - vx(1,0) = vx(1,0) + Nwx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nwx(1,a)*yl(j,a); - vx(1,2) = vx(1,2) + Nwx(2,a)*yl(j,a); + // Active stress activation along fiber, sheet and sheet-normal + const double ya_g_f = eigen_view(ya_l_f).dot(Nwm); + const double ya_g_s = eigen_view(ya_l_s).dot(Nwm); + const double ya_g_n = eigen_view(ya_l_n).dot(Nwm); - vx(2,0) = vx(2,0) + Nwx(0,a)*yl(k,a); - vx(2,1) = vx(2,1) + Nwx(1,a)*yl(k,a); - vx(2,2) = vx(2,2) + Nwx(2,a)*yl(k,a); + // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) + const Matrix<3> vx = vel * Nwxm.transpose(); + const Matrix<3> F = Matrix<3>::Identity() + disp * Nwxm.transpose(); - F(0,0) = F(0,0) + Nwx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nwx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nwx(2,a)*dl(i,a); - - F(1,0) = F(1,0) + Nwx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nwx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nwx(2,a)*dl(j,a); - - F(2,0) = F(2,0) + Nwx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nwx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nwx(2,a)*dl(k,a); - - ya_g_f = ya_g_f + Nw(a) * ya_l_f(a); - ya_g_s = ya_g_s + Nw(a) * ya_l_s(a); - ya_g_n = ya_g_n + Nw(a) * ya_l_n(a); - } - - double Jac = mat_fun::mat_det(F, 3); - auto Fi = mat_fun::mat_inv(F, 3); + double Jac = F.determinant(); + const Matrix<3> Fi = F.inverse(); // Pressure and its time derivative // @@ -1283,18 +1108,17 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Compute deviatoric 2nd Piola-Kirchhoff stress tensor (Siso) and // isochoric elasticity tensor in Voigt notation (Dm) // - Array Siso(3,3), Dm(6,6); + Matrix<3> Siso; + Matrix<6> Dm; double Ja = 0; - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g_f, + mat_models::compute_pk2cc<3>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<3>(fN), ya_g_f, ya_g_s, ya_g_n, Siso, Dm, Ja); - // Viscous 2nd Piola-Kirchhoff stress and tangent contributions - Array Svis(3,3); - Array3 Kvis_u(9, eNoNw, eNoNw); - Array3 Kvis_v(9, eNoNw, eNoNw); - - mat_models::compute_visc_stress_and_tangent(dmn, eNoNw, Nwx, vx, F, Svis, Kvis_u, Kvis_v); - + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. Reuse the + // previous Gauss point's when shape function gradients are constant within an + // element (e.g. linear triangles and tetrahedra). + static mat_models::ViscousResponse<3> visc; + visc.update(dmn, eNoNw, Nwx, vx, F, recompute_visc); // Compute rho and beta depending on the volumetric penalty model // @@ -1316,260 +1140,195 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, } // Total isochoric 2nd Piola-Kirchhoff stress (Elastic + Viscous) - Siso = Siso + Svis; + Siso += visc.S(); // Deviatoric 1st Piola-Kirchhoff tensor (P) // - auto Pdev = mat_fun::mat_mul(F, Siso); + const Matrix<3> Pdev = F * Siso; // Shape function gradients in the current configuration // - Array NxFi(3,eNoNw); - - for (int a = 0; a < eNoNw; a++) { - NxFi(0,a) = Nwx(0,a)*Fi(0,0) + Nwx(1,a)*Fi(1,0) + Nwx(2,a)*Fi(2,0); - NxFi(1,a) = Nwx(0,a)*Fi(0,1) + Nwx(1,a)*Fi(1,1) + Nwx(2,a)*Fi(2,1); - NxFi(2,a) = Nwx(0,a)*Fi(0,2) + Nwx(1,a)*Fi(1,2) + Nwx(2,a)*Fi(2,2); - } + const NodalMatrix<3> NxFi = Fi.transpose() * Nwxm; // Velocity gradient in current configuration - auto VxFi = mat_mul(vx, Fi); + const Matrix<3> VxFi = vx * Fi; double rC = beta*pd + VxFi(0,0) + VxFi(1,1) + VxFi(2,2); double rCl = -p + tauC*rC; - // Local residual - // - double T1, T2, T3; - - for (int a = 0; a < eNoNw; a++) { - T1 = Jac*rho*vd(0)*Nw(a); - T2 = Pdev(0,0)*Nwx(0,a) + Pdev(0,1)*Nwx(1,a) + Pdev(0,2)*Nwx(2,a); - T3 = Jac*rCl*NxFi(0,a); - lR(0,a) = lR(0,a) + w*(T1 + T2 + T3); - - T1 = Jac*rho*vd(1)*Nw(a); - T2 = Pdev(1,0)*Nwx(0,a) + Pdev(1,1)*Nwx(1,a) + Pdev(1,2)*Nwx(2,a); - T3 = Jac*rCl*NxFi(1,a); - lR(1,a) = lR(1,a) + w*(T1 + T2 + T3); - - T1 = Jac*rho*vd(2)*Nw(a); - T2 = Pdev(2,0)*Nwx(0,a) + Pdev(2,1)*Nwx(1,a) + Pdev(2,2)*Nwx(2,a); - T3 = Jac*rCl*NxFi(2,a); - lR(2,a) = lR(2,a) + w*(T1 + T2 + T3); - } + // Inertia, the divergence of Pdev, and the pressure/volumetric term + lRv += w * (Jac*rho * vd * Nwm.transpose() + Pdev * Nwxm + Jac*rCl * NxFi); - // Auxilary quantities for computing stiffness tensors + // Strain-displacement matrix; Bm[a] maps node a to Voigt strain // - Array3 Bm(6,3,eNoNw); + std::array, consts::maxNoN> Bm; - for (int a = 0; a < eNoNw; a++) { - Bm(0,0,a) = Nwx(0,a)*F(0,0); - Bm(0,1,a) = Nwx(0,a)*F(1,0); - Bm(0,2,a) = Nwx(0,a)*F(2,0); - - Bm(1,0,a) = Nwx(1,a)*F(0,1); - Bm(1,1,a) = Nwx(1,a)*F(1,1); - Bm(1,2,a) = Nwx(1,a)*F(2,1); - - Bm(2,0,a) = Nwx(2,a)*F(0,2); - Bm(2,1,a) = Nwx(2,a)*F(1,2); - Bm(2,2,a) = Nwx(2,a)*F(2,2); - - Bm(3,0,a) = (Nwx(0,a)*F(0,1) + F(0,0)*Nwx(1,a)); - Bm(3,1,a) = (Nwx(0,a)*F(1,1) + F(1,0)*Nwx(1,a)); - Bm(3,2,a) = (Nwx(0,a)*F(2,1) + F(2,0)*Nwx(1,a)); - - Bm(4,0,a) = (Nwx(1,a)*F(0,2) + F(0,1)*Nwx(2,a)); - Bm(4,1,a) = (Nwx(1,a)*F(1,2) + F(1,1)*Nwx(2,a)); - Bm(4,2,a) = (Nwx(1,a)*F(2,2) + F(2,1)*Nwx(2,a)); - - Bm(5,0,a) = (Nwx(2,a)*F(0,0) + F(0,2)*Nwx(0,a)); - Bm(5,1,a) = (Nwx(2,a)*F(1,0) + F(1,2)*Nwx(0,a)); - Bm(5,2,a) = (Nwx(2,a)*F(2,0) + F(2,2)*Nwx(0,a)); - } - - Array VxNx(3,eNoNw); + const Matrix<3> Ft = F.transpose(); for (int a = 0; a < eNoNw; a++) { - VxNx(0,a) = VxFi(0,0)*NxFi(0,a) + VxFi(1,0)*NxFi(1,a) + VxFi(2,0)*NxFi(2,a); - VxNx(1,a) = VxFi(0,1)*NxFi(0,a) + VxFi(1,1)*NxFi(1,a) + VxFi(2,1)*NxFi(2,a); - VxNx(2,a) = VxFi(0,2)*NxFi(0,a) + VxFi(1,2)*NxFi(1,a) + VxFi(2,2)*NxFi(2,a); + const auto g = Nwxm.col(a); // grad(N_a) + + Bm[a].row(0) = g(0) * Ft.row(0); // dE_11 + Bm[a].row(1) = g(1) * Ft.row(1); // dE_22 + Bm[a].row(2) = g(2) * Ft.row(2); // dE_33 + Bm[a].row(3) = g(0) * Ft.row(1) + g(1) * Ft.row(0); // 2 dE_12 + Bm[a].row(4) = g(1) * Ft.row(2) + g(2) * Ft.row(1); // 2 dE_23 + Bm[a].row(5) = g(2) * Ft.row(0) + g(0) * Ft.row(2); // 2 dE_31 } + const NodalMatrix<3> VxNx = VxFi.transpose() * NxFi; + // Tangent (stiffness) matrices // double r13 = 1.0 / 3.0; double r23 = 2.0 / 3.0; - double NxSNx{0.0}, BtDB{0.0}; + double NxSNx{0.0}, BtDB{0.0}, T1{0.0}, T2{0.0}, T3{0.0}; double Tv{0.0}, Ku{0.0}; - Array DBm(6,3); - for (int b = 0; b < eNoNw; b++) { - mat_mul(Dm, Bm.rslice(b), DBm); + const Eigen::Matrix DBm = Dm * Bm[b]; + + // Geometric stiffness: Siso*grad(N_b) + const Eigen::Vector3d SisoNx = Siso * Nwxm.col(b); for (int a = 0; a < eNoNw; a++) { - NxSNx = Nwx(0,a)*Siso(0,0)*Nwx(0,b) - + Nwx(0,a)*Siso(0,1)*Nwx(1,b) + Nwx(0,a)*Siso(0,2)*Nwx(2,b) - + Nwx(1,a)*Siso(1,0)*Nwx(0,b) + Nwx(1,a)*Siso(1,1)*Nwx(1,b) - + Nwx(1,a)*Siso(1,2)*Nwx(2,b) + Nwx(2,a)*Siso(2,0)*Nwx(0,b) - + Nwx(2,a)*Siso(2,1)*Nwx(1,b) + Nwx(2,a)*Siso(2,2)*Nwx(2,b); + NxSNx = Nwxm.col(a).dot(SisoNx); // dM1_dV1 + af/am *dM_1/dU_1 - BtDB = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + - Bm(2,0,a)*DBm(2,0) + Bm(3,0,a)*DBm(3,0) + - Bm(4,0,a)*DBm(4,0) + Bm(5,0,a)*DBm(5,0); + BtDB = Bm[a].col(0).dot(DBm.col(0)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(0,b); - Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(0,a,b)); - lKd(0,a,b) = lKd(0,a,b) + Ku; + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(0,a,b)); + lKd(0,a,b) += Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(0,a)*NxFi(0,b); - Tv = af*Kvis_v(0,a,b); - lK(0,a,b) = lK(0,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(0,a,b); + lK(0,a,b) += w*(T2 + Tv) + afm*Ku; // dM_1/dV_2 + af/am *dM_1/dU_2 - BtDB = Bm(0,0,a)*DBm(0,1) + Bm(1,0,a)*DBm(1,1) + - Bm(2,0,a)*DBm(2,1) + Bm(3,0,a)*DBm(3,1) + - Bm(4,0,a)*DBm(4,1) + Bm(5,0,a)*DBm(5,1); + BtDB = Bm[a].col(0).dot(DBm.col(1)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(1,b) - NxFi(1,a)*NxFi(0,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(1,a,b)); - lKd(1,a,b) = lKd(1,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(1,a,b)); + lKd(1,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(1,b); - Tv = af*Kvis_v(1,a,b); - lK(1,a,b) = lK(1,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(1,a,b); + lK(1,a,b) += w*(T2 + Tv) + afm*Ku; // dM_1/dV_3 + af/am *dM_1/dU_3 // - BtDB = Bm(0,0,a)*DBm(0,2) + Bm(1,0,a)*DBm(1,2) + - Bm(2,0,a)*DBm(2,2) + Bm(3,0,a)*DBm(3,2) + - Bm(4,0,a)*DBm(4,2) + Bm(5,0,a)*DBm(5,2); + BtDB = Bm[a].col(0).dot(DBm.col(2)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(2,b) - NxFi(2,a)*NxFi(0,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(2,a,b)); - lKd(2,a,b) = lKd(2,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(2,a,b)); + lKd(2,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(2,b); - Tv = af*Kvis_v(2,a,b); - lK(2,a,b) = lK(2,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(2,a,b); + lK(2,a,b) += w*(T2 + Tv) + afm*Ku; // dM_2/dV_1 + af/am *dM_2/dU_1 // - BtDB = Bm(0,1,a)*DBm(0,0) + Bm(1,1,a)*DBm(1,0) + - Bm(2,1,a)*DBm(2,0) + Bm(3,1,a)*DBm(3,0) + - Bm(4,1,a)*DBm(4,0) + Bm(5,1,a)*DBm(5,0); + BtDB = Bm[a].col(1).dot(DBm.col(0)); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(0,b) - NxFi(0,a)*NxFi(1,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(3,a,b)); - lKd(3,a,b) = lKd(3,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(3,a,b)); + lKd(3,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(0,b); - Tv = af*Kvis_v(3,a,b); + Tv = af*visc.dv(3,a,b); - lK(4,a,b) = lK(4,a,b) + w*(T2 + Tv) + afm*Ku; + lK(4,a,b) += w*(T2 + Tv) + afm*Ku; // dM_2/dV_2 + af/am *dM_2/dU_2 // - BtDB = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + - Bm(2,1,a)*DBm(2,1) + Bm(3,1,a)*DBm(3,1) + - Bm(4,1,a)*DBm(4,1) + Bm(5,1,a)*DBm(5,1); + BtDB = Bm[a].col(1).dot(DBm.col(1)); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(1,b); - Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(4,a,b)); - lKd(4,a,b) = lKd(4,a,b) + Ku; + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(4,a,b)); + lKd(4,a,b) += Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(1,a)*NxFi(1,b); - Tv = af*Kvis_v(4,a,b); - lK(5,a,b) = lK(5,a,b) + w*(T2 + Tv) + afm*Ku; - + Tv = af*visc.dv(4,a,b); + lK(5,a,b) += w*(T2 + Tv) + afm*Ku; // dM_2/dV_3 + af/am *dM_2/dU_3 // - BtDB = Bm(0,1,a)*DBm(0,2) + Bm(1,1,a)*DBm(1,2) + - Bm(2,1,a)*DBm(2,2) + Bm(3,1,a)*DBm(3,2) + - Bm(4,1,a)*DBm(4,2) + Bm(5,1,a)*DBm(5,2); + BtDB = Bm[a].col(1).dot(DBm.col(2)); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(2,b) - NxFi(2,a)*NxFi(1,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(5,a,b)); - lKd(5,a,b) = lKd(5,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(5,a,b)); + lKd(5,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(2,b); - Tv = af*Kvis_v(5,a,b); - lK(6,a,b) = lK(6,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(5,a,b); + lK(6,a,b) += w*(T2 + Tv) + afm*Ku; // dM_3/dV_1 + af/am *dM_3/dU_1 // - BtDB = Bm(0,2,a)*DBm(0,0) + Bm(1,2,a)*DBm(1,0) + - Bm(2,2,a)*DBm(2,0) + Bm(3,2,a)*DBm(3,0) + - Bm(4,2,a)*DBm(4,0) + Bm(5,2,a)*DBm(5,0); + BtDB = Bm[a].col(2).dot(DBm.col(0)); T1 = Jac*rho*vd(2)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(0,b) - NxFi(0,a)*NxFi(2,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(6,a,b)); - lKd(6,a,b) = lKd(6,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(6,a,b)); + lKd(6,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(0,b); - Tv = af*Kvis_v(6,a,b); - lK(8,a,b) = lK(8,a,b) + w*(T2 + Tv) + afm*Ku; + Tv = af*visc.dv(6,a,b); + lK(8,a,b) += w*(T2 + Tv) + afm*Ku; // dM_3/dV_2 + af/am *dM_3/dU_2 // - BtDB = Bm(0,2,a)*DBm(0,1) + Bm(1,2,a)*DBm(1,1) + - Bm(2,2,a)*DBm(2,1) + Bm(3,2,a)*DBm(3,1) + - Bm(4,2,a)*DBm(4,1) + Bm(5,2,a)*DBm(5,1); + BtDB = Bm[a].col(2).dot(DBm.col(1)); T1 = Jac*rho*vd(2)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(1,b) - NxFi(1,a)*NxFi(2,b)); - Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(7,a,b)); - lKd(7,a,b) = lKd(7,a,b) + Ku; + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(7,a,b)); + lKd(7,a,b) += Ku; T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(1,b); - Tv = af*Kvis_v(7,a,b); + Tv = af*visc.dv(7,a,b); - lK(9,a,b) = lK(9,a,b) + w*(T2 + Tv) + afm*Ku; + lK(9,a,b) += w*(T2 + Tv) + afm*Ku; // dM_3/dV_3 + af/am *dM_3/dU_3 // - BtDB = Bm(0,2,a)*DBm(0,2) + Bm(1,2,a)*DBm(1,2) + - Bm(2,2,a)*DBm(2,2) + Bm(3,2,a)*DBm(3,2) + - Bm(4,2,a)*DBm(4,2) + Bm(5,2,a)*DBm(5,2); + BtDB = Bm[a].col(2).dot(DBm.col(2)); T1 = Jac*rho*vd(2)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(2,b); - Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(8,a,b)); - lKd(8,a,b) = lKd(8,a,b) + Ku; + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(8,a,b)); + lKd(8,a,b) += Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(2,a)*NxFi(2,b); - Tv = af*Kvis_v(8,a,b); + Tv = af*visc.dv(8,a,b); - lK(10,a,b) = lK(10,a,b) + w*(T2 + Tv) + afm*Ku; + lK(10,a,b) += w*(T2 + Tv) + afm*Ku; } } @@ -1580,19 +1339,21 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // dM_0/dP T0 = am*tauC*beta + af*(tauC*dbeta*pd - 1.0); T1 = T0*NxFi(0,a)*Nq(b) + af*drho*vd(0)*Nw(a)*Nq(b); - lK(3,a,b) = lK(3,a,b) + w*Jac*T1; + lK(3,a,b) += w*Jac*T1; // dM_1/dP T1 = T0*NxFi(1,a)*Nq(b) + af*drho*vd(1)*Nw(a)*Nq(b); - lK(7,a,b) = lK(7,a,b) + w*Jac*T1; + lK(7,a,b) += w*Jac*T1; // dM_2/dP T1 = T0*NxFi(2,a)*Nq(b) + af*drho*vd(2)*Nw(a)*Nq(b); - lK(11,a,b) = lK(11,a,b) + w*Jac*T1; + lK(11,a,b) += w*Jac*T1; } } } +/// @brief Replicates 'SUBROUTINE USTRUCT_DOASSEM(d, eqN, lKd, lK, lR)' +// /// @brief Replicates 'SUBROUTINE USTRUCT_DOASSEM(d, eqN, lKd, lK, lR)' // void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, diff --git a/Code/Source/solver/ustruct.h b/Code/Source/solver/ustruct.h index 953d4d035..c89e99c61 100644 --- a/Code/Source/solver/ustruct.h +++ b/Code/Source/solver/ustruct.h @@ -37,7 +37,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &dl, const Array &bfl, const Array &fN, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK, Array3 &lKd); + Array &lR, Array3 &lK, Array3 &lKd, + const bool recompute_visc); void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, @@ -53,7 +54,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &dl, const Array &bfl, const Array &fN, const Vector &ya_l_f, const Vector &ya_l_s, const Vector &ya_l_n, - Array &lR, Array3 &lK, Array3 &lKd); + Array &lR, Array3 &lK, Array3 &lKd, + const bool recompute_visc); void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, const Array3& lK, const Array& lR);