diff --git a/Code/Source/solver/mat_fun.cpp b/Code/Source/solver/mat_fun.cpp index 582db47da..bb0a4c081 100644 --- a/Code/Source/solver/mat_fun.cpp +++ b/Code/Source/solver/mat_fun.cpp @@ -507,15 +507,10 @@ 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 + // Fall back onto the overload of this function that takes the result matrix + // as an output argument. + mat_mul(A, B, result); return result; } @@ -524,8 +519,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 +949,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..caef69cea 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -1692,10 +1692,11 @@ void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array< 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); } + mat_mul(ddev, Nx_Fi, ddev_Nx_Fi); + mat_mul(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); 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) +