diff --git a/Code/Source/solver/mat_fun.cpp b/Code/Source/solver/mat_fun.cpp index 582db47da..d2e4c6592 100644 --- a/Code/Source/solver/mat_fun.cpp +++ b/Code/Source/solver/mat_fun.cpp @@ -507,15 +507,9 @@ mat_mul(const Array& A, const Array& B) Array result(A_num_rows, B_num_cols); - for (int i = 0; i < A_num_rows; i++) { - for (int j = 0; j < B_num_cols; j++) { - double sum = 0.0; - for (int k = 0; k < A_num_cols; k++) { - sum += A(i,k) * B(k,j); - } - result(i,j) = sum; - } - } + // Delegate rather than repeat the loop, so callers of this form reach the + // fixed-shape fast paths too. Most call sites use this overload. + mat_mul(A, B, result); return result; } @@ -524,8 +518,62 @@ mat_mul(const Array& A, const Array& B) /// /// Compute result directly into the passed argument. // +namespace { + +/// @brief Fixed-shape matrix product, C = A*B, mapped onto the existing buffers. +/// +/// Used when sizes are known at compile time. +template +inline void mat_mul_fixed(const Array& A, const Array& B, Array& C) +{ + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data()); + Eigen::Map> c(C.data()); + + c.noalias() = a * b; +} + +/// @brief As mat_mul_fixed, but with the column count known only at run time. +/// +/// Used where the right operand has one column per element node, so its width +/// depends on the element type. The row counts are still compile-time, which is +/// where most of the benefit comes from. +template +inline void mat_mul_fixed_rows(const Array& A, const Array& B, Array& C) +{ + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data(), K, B.ncols()); + Eigen::Map> c(C.data(), M, C.ncols()); + + c.noalias() = a * b; +} + +} // namespace + void mat_mul(const Array& A, const Array& B, Array& result) { + // Fixed-shape fast paths for the products that dominate the element loops. + // + // 3x3 * 3x3 F*S in struct_3d; vx*Fi, ddev*Fit, Fi*ddev_Fit and the + // potential-viscosity products in mat_models; F^T*F in cep + // 3x3 * 3xeNoN ddev*Nx_Fi and vx_Fi*Nx_Fi in the viscous tangent + // 6x6 * 6x3 the material stiffness product D*B in struct_3d/ustruct_3d_m + if (A.nrows() == 3 && A.ncols() == 3 && B.nrows() == 3 && + result.nrows() == 3 && result.ncols() == B.ncols()) { + if (B.ncols() == 3) { + mat_mul_fixed<3, 3, 3>(A, B, result); + } else { + mat_mul_fixed_rows<3, 3>(A, B, result); + } + return; + } + + if (A.nrows() == 6 && A.ncols() == 6 && B.nrows() == 6 && B.ncols() == 3 && + result.nrows() == 6 && result.ncols() == 3) { + mat_mul_fixed<6, 6, 3>(A, B, result); + return; + } + int A_num_rows = A.nrows(); int A_num_cols = A.ncols(); int B_num_rows = B.nrows(); @@ -900,36 +948,6 @@ transpose(const Array& A) return result; } -void mat_mul6x3(const Array& A, const Array& B, Array& C) -{ - #define mat_mul6x3_unroll - #ifdef mat_mul6x3_unroll - auto a = A.data(); - auto b = B.data(); - auto c = C.data(); - - c[0] = a[0]*b[0] + a[6]*b[1] + a[12]*b[2] + a[18]*b[3] + a[24]*b[4] + a[30]*b[5]; - c[1] = a[1]*b[0] + a[7]*b[1] + a[13]*b[2] + a[19]*b[3] + a[25]*b[4] + a[31]*b[5]; - c[2] = a[2]*b[0] + a[8]*b[1] + a[14]*b[2] + a[20]*b[3] + a[26]*b[4] + a[32]*b[5]; - c[3] = a[3]*b[0] + a[9]*b[1] + a[15]*b[2] + a[21]*b[3] + a[27]*b[4] + a[33]*b[5]; - c[4] = a[4]*b[0] + a[10]*b[1] + a[16]*b[2] + a[22]*b[3] + a[28]*b[4] + a[34]*b[5]; - c[5] = a[5]*b[0] + a[11]*b[1] + a[17]*b[2] + a[23]*b[3] + a[29]*b[4] + a[35]*b[5]; - - c[6] = a[0]*b[6] + a[6]*b[7] + a[12]*b[8] + a[18]*b[9] + a[24]*b[10] + a[30]*b[11]; - c[7] = a[1]*b[6] + a[7]*b[7] + a[13]*b[8] + a[19]*b[9] + a[25]*b[10] + a[31]*b[11]; - c[8] = a[2]*b[6] + a[8]*b[7] + a[14]*b[8] + a[20]*b[9] + a[26]*b[10] + a[32]*b[11]; - c[9] = a[3]*b[6] + a[9]*b[7] + a[15]*b[8] + a[21]*b[9] + a[27]*b[10] + a[33]*b[11]; - c[10] = a[4]*b[6] + a[10]*b[7] + a[16]*b[8] + a[22]*b[9] + a[28]*b[10] + a[34]*b[11]; - c[11] = a[5]*b[6] + a[11]*b[7] + a[17]*b[8] + a[23]*b[9] + a[29]*b[10] + a[35]*b[11]; - - c[12] = a[0]*b[12] + a[6]*b[13] + a[12]*b[14] + a[18]*b[15] + a[24]*b[16] + a[30]*b[17]; - c[13] = a[1]*b[12] + a[7]*b[13] + a[13]*b[14] + a[19]*b[15] + a[25]*b[16] + a[31]*b[17]; - c[14] = a[2]*b[12] + a[8]*b[13] + a[14]*b[14] + a[20]*b[15] + a[26]*b[16] + a[32]*b[17]; - c[15] = a[3]*b[12] + a[9]*b[13] + a[15]*b[14] + a[21]*b[15] + a[27]*b[16] + a[33]*b[17]; - c[16] = a[4]*b[12] + a[10]*b[13] + a[16]*b[14] + a[22]*b[15] + a[28]*b[16] + a[34]*b[17]; - c[17] = a[5]*b[12] + a[11]*b[13] + a[17]*b[14] + a[23]*b[15] + a[29]*b[16] + a[35]*b[17]; - #endif -} }; diff --git a/Code/Source/solver/mat_fun.h b/Code/Source/solver/mat_fun.h index 5428479c7..caf83b628 100644 --- a/Code/Source/solver/mat_fun.h +++ b/Code/Source/solver/mat_fun.h @@ -97,7 +97,6 @@ namespace mat_fun { Vector mat_mul(const Array& A, const Vector& v); Array mat_mul(const Array& A, const Array& B); void mat_mul(const Array& A, const Array& B, Array& result); - void mat_mul6x3(const Array& A, const Array& B, Array& C); Array mat_symm(const Array& A, const int nd); Array mat_symm_prod(const Vector& u, const Vector& v, const int nd); diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index cda281a68..73133f518 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -1559,6 +1559,117 @@ void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, } } +namespace { + +// --------------------------------------------------------------------------- +// Fixed-size stand-ins for the mat_fun helpers the viscous kernels use. +// +// Only the storage changes -- from heap-allocated Array to stack-resident +// Eigen -- so the arithmetic below is transcribed operation-for-operation from +// mat_fun.cpp and stays bit-identical. Two of the originals are worth calling +// out: mat_det() expands a 3x3 determinant recursively, allocating a minor and +// calling pow() three times per evaluation, and mat_dev() costs four +// allocations because Array's operator- and operator* both return by value. +// +// eNoN depends on element type; HEX27 is the largest, so the max-size Eigen +// form keeps the nsd-by-eNoN scratch on the stack while its column count stays +// a run-time value. +// --------------------------------------------------------------------------- + +constexpr int MAX_ENON = 27; + +template +struct VType { + using Mat = Eigen::Matrix; + using MatX = Eigen::Matrix; +}; + +inline double det_fixed(const Eigen::Matrix& A) +{ + return A(0,0)*A(1,1) - A(0,1)*A(1,0); +} + +/// @brief mat_det()'s cofactor expansion along row 0, accumulated in its order. +/// +/// The original multiplies each term by pow(-1.0, 2+i), which is exactly +1, +/// -1, +1, so the signs are folded in here without changing any rounding. +inline double det_fixed(const Eigen::Matrix& A) +{ + double D = 0.0; + D = D + A(0,0) * (A(1,1)*A(2,2) - A(1,2)*A(2,1)); + D = D - A(0,1) * (A(1,0)*A(2,2) - A(1,2)*A(2,0)); + D = D + A(0,2) * (A(1,0)*A(2,1) - A(1,1)*A(2,0)); + return D; +} + +/// @brief mat_inv() for the fixed sizes, dividing each cofactor by the +/// determinant exactly as the original does. +/// +/// Eigen's .inverse() computes 1/d once and multiplies through, which rounds +/// differently, so it deliberately is not used here. +inline Eigen::Matrix inv_fixed(const Eigen::Matrix& A, + const double d) +{ + Eigen::Matrix Ai; + Ai(0,0) = A(1,1) / d; + Ai(0,1) = -A(0,1) / d; + Ai(1,0) = -A(1,0) / d; + Ai(1,1) = A(0,0) / d; + return Ai; +} + +inline Eigen::Matrix inv_fixed(const Eigen::Matrix& A, + const double d) +{ + Eigen::Matrix Ai; + Ai(0,0) = (A(1,1)*A(2,2) - A(1,2)*A(2,1)) / d; + Ai(0,1) = (A(0,2)*A(2,1) - A(0,1)*A(2,2)) / d; + Ai(0,2) = (A(0,1)*A(1,2) - A(0,2)*A(1,1)) / d; + + Ai(1,0) = (A(1,2)*A(2,0) - A(1,0)*A(2,2)) / d; + Ai(1,1) = (A(0,0)*A(2,2) - A(0,2)*A(2,0)) / d; + Ai(1,2) = (A(0,2)*A(1,0) - A(0,0)*A(1,2)) / d; + + Ai(2,0) = (A(1,0)*A(2,1) - A(1,1)*A(2,0)) / d; + Ai(2,1) = (A(0,1)*A(2,0) - A(0,0)*A(2,1)) / d; + Ai(2,2) = (A(0,0)*A(1,1) - A(0,1)*A(1,0)) / d; + return Ai; +} + +/// @brief mat_symm(): S(i,j) = 0.5 * (A(i,j) + A(j,i)). +template +inline Eigen::Matrix symm_fixed(const Eigen::Matrix& A) +{ + Eigen::Matrix S; + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + S(i,j) = 0.5 * (A(i,j) + A(j,i)); + } + } + return S; +} + +/// @brief mat_dev(): A - (tr(A) / nsd) * I, elementwise as the original. +template +inline Eigen::Matrix dev_fixed(const Eigen::Matrix& A) +{ + double trA = 0.0; + for (int i = 0; i < nsd; i++) { + trA = trA + A(i,i); + } + const double c = trA / static_cast(nsd); + + Eigen::Matrix R; + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + R(i,j) = A(i,j) - c * (i == j ? 1.0 : 0.0); + } + } + return R; +} + +} // namespace + /** * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid * with a viscous pseudo-potential model. @@ -1580,44 +1691,41 @@ void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, * @param Kvis_u Viscous tangent matrix contribution due to displacement * @param Kvis_v Visous 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; +namespace { - // Number of spatial dimensions - int nsd = F.nrows(); +/// @brief compute_visc_stress_potential() with the spatial dimension fixed at +/// compile time. Same transformation as visc_newtonian_impl: heap temporaries +/// become stack ones, arithmetic unchanged. +template +void visc_potential_impl(const double mu, const int eNoN, const Array& Nx, + const Array& vx, const Array& F, + Array& Svis, Array3& Kvis_u, + Array3& Kvis_v) +{ + using Mat = typename VType::Mat; + using MatX = typename VType::MatX; // Initialize Svis, Kvis_u, Kvis_v to zero Svis = 0.0; Kvis_u = 0.0; Kvis_v = 0.0; + Eigen::Map Fm(F.data()); + Eigen::Map vxm(vx.data()); + Eigen::Map> Nxm(Nx.data(), nsd, eNoN); // 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 Mat F_Ft = Fm * Fm.transpose(); + const Mat Ft_vx = Fm.transpose() * vxm; + const Mat F_vxt = Fm * vxm.transpose(); + + // F_Nx(i,a) = sum_j F(i,j) * Nx(j,a); vx_Nx likewise. + const MatX F_Nx = Fm * Nxm; + const MatX vx_Nx = vxm * 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); + Eigen::Map(Svis.data()) = mu * symm_fixed(Ft_vx); // Tangent matrix contributions due to viscosity for (int b = 0; b < eNoN; ++b) { @@ -1638,6 +1746,18 @@ void compute_visc_stress_potential(const double mu, const int eNoN, const Array< } } +} // namespace + +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) { + if (F.nrows() == 3) { + visc_potential_impl<3>(mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + } else { + visc_potential_impl<2>(mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + } +} + + /** * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid * with a Newtonian fluid-like viscosity model. @@ -1659,49 +1779,56 @@ void compute_visc_stress_potential(const double mu, const int eNoN, const Array< * @param Kvis_u Viscous tangent matrix contribution due to displacement * @param Kvis_v Visous 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; +namespace { - // Number of spatial dimensions - int nsd = F.nrows(); +/// @brief compute_visc_stress_newtonian() with the spatial dimension fixed at +/// compile time, so every temporary lives on the stack instead of the heap. +/// +/// The original allocated roughly fifteen Arrays per Gauss point: one for each +/// named intermediate, plus four hidden inside mat_dev() and one inside every +/// mat_det() call. All of that is gone; the arithmetic is unchanged. +template +void visc_newtonian_impl(const double mu, const int eNoN, const Array& Nx, + const Array& vx, const Array& F, + Array& Svis, Array3& Kvis_u, + Array3& Kvis_v) +{ + using Mat = typename VType::Mat; + using MatX = typename VType::MatX; // Initialize Svis, Kvis_u, Kvis_v to zero Svis = 0.0; Kvis_u = 0.0; Kvis_v = 0.0; - // Get identity matrix, Jacobian, and F^-1 - auto Idm = mat_id(nsd); - auto J = mat_det(F, nsd); - auto Fi = mat_inv(F, nsd); + // Array and Eigen are both column-major, so these alias the caller's + // buffers rather than copying them. + Eigen::Map Fm(F.data()); + Eigen::Map vxm(vx.data()); + Eigen::Map> Nxm(Nx.data(), nsd, eNoN); + + // Jacobian and F^-1 + const double J = det_fixed(Mat(Fm)); + if (utils::is_zero(fabs(J))) { + throw std::runtime_error("Singular matrix detected to compute inverse"); + } + const Mat Fi = inv_fixed(Mat(Fm), J); // 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); + // vx_Fi: Velocity gradient in current configuration + const Mat vx_Fi = vxm * Fi; // 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); - } - } - ddev_Nx_Fi = mat_mul(ddev, Nx_Fi); - vx_Fi_Nx_Fi = mat_mul(vx_Fi, Nx_Fi); - } + const Mat ddev = dev_fixed(symm_fixed(vx_Fi)); + + // Nx_Fi(i,a) = sum_j Nx(j,a) * Fi(j,i), i.e. Fi^T * Nx. + const MatX Nx_Fi = Fi.transpose() * Nxm; + const MatX ddev_Nx_Fi = ddev * Nx_Fi; + const MatX 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; + const Mat Fi_ddev_Fit = Fi * (ddev * Fi.transpose()); + Eigen::Map(Svis.data()) = 2.0 * mu * J * Fi_ddev_Fit; // Tangent matrix contributions due to viscosity double r2d = 2.0 / nsd; @@ -1717,13 +1844,14 @@ void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array< 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))); // Derivative of the residual w.r.t velocity - Kvis_v(ii,a,b) = mu * J * (Nx_Fi_Nx_Fi * Idm(i,j) + + // Idm(i,j) was a freshly allocated identity matrix per call. + Kvis_v(ii,a,b) = mu * J * (Nx_Fi_Nx_Fi * (i == j ? 1.0 : 0.0) + Nx_Fi(i,b) * Nx_Fi(j,a) - r2d * Nx_Fi(i,a) * Nx_Fi(j,b)); } } @@ -1731,6 +1859,18 @@ void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array< } } +} // namespace + +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) { + // Dispatch on the spatial dimension; both 2D and 3D solid kernels call here. + if (F.nrows() == 3) { + visc_newtonian_impl<3>(mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + } else { + visc_newtonian_impl<2>(mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + } +} + /** * @brief Get the solid viscous PK2 stress and corresponding tangent matrix contributions diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index 9ab5feca6..41848e485 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -735,6 +735,9 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, for (int b = 0; b < eNoN; b++) { + // Material stiffness (D*B) + mat_mul(Dm, Bm.rslice(b), DBm); + for (int a = 0; a < eNoN; a++) { // Geometric stiffness @@ -746,9 +749,6 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, T1 = amd*N(a)*N(b) + afu*NxSNx; - // Material Stiffness (Bt*D*B) - mat_mul(Dm, Bm.rslice(b), DBm); - // dM1/du1 // Material stiffness: Bt*D*B BmDBm = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 69be045c7..a698050e7 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -1403,7 +1403,12 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, double NxSNx{0.0}, BtDB{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); + 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) @@ -1411,8 +1416,6 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + 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); - auto DBm = mat_mul(Dm, Bm.rslice(b)); - // 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) +