diff --git a/Code/Source/svFSI/fluid.cpp b/Code/Source/svFSI/fluid.cpp index 2acdcec..530db7e 100644 --- a/Code/Source/svFSI/fluid.cpp +++ b/Code/Source/svFSI/fluid.cpp @@ -577,7 +577,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag for (int g = 0; g < fs[0].nG; g++) { if (g == 0 || !fs[1].lShpF) { - auto Nx = fs[1].Nx.slice(g); + auto Nx = fs[1].Nx.rslice(g); nn::gnn(fs[1].eNoN, nsd, nsd, Nx, xql, Nqx, Jac, ksix); if (utils::is_zero(Jac)) { throw std::runtime_error("[construct_fluid] Jacobian for element " + std::to_string(e) + " is < 0."); @@ -585,15 +585,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag } if (g == 0 || !fs[0].lShpF) { - auto Nx = fs[0].Nx.slice(g); + auto Nx = fs[0].Nx.rslice(g); nn::gnn(fs[0].eNoN, nsd, nsd, Nx, xwl, Nwx, Jac, ksix); if (utils::is_zero(Jac)) { throw std::runtime_error("[construct_fluid] Jacobian for element " + std::to_string(e) + " is < 0."); } if (!vmsStab) { - auto Nx = fs[0].Nx.slice(g); - auto Nxx = fs[0].Nxx.slice(g); + auto Nx = fs[0].Nx.rslice(g); + auto Nxx = fs[0].Nxx.rslice(g); nn::gn_nxx(l, fs[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx); } } @@ -603,14 +603,16 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // Compute momentum residual and tangent matrix. // if (nsd == 3) { - auto N0 = fs[0].N.col(g); - auto N1 = fs[1].N.col(g); - fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); + auto N0 = fs[0].N.rcol(g); + auto N1 = fs[1].N.rcol(g); + fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); } else if (nsd == 2) { - auto N0 = fs[0].N.col(g); - auto N1 = fs[1].N.col(g); - fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); + auto N0 = fs[0].N.rcol(g); + auto N1 = fs[1].N.rcol(g); + fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); } } // g: loop @@ -622,7 +624,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // for (int g = 0; g < fs[1].nG; g++) { if (g == 0 || !fs[0].lShpF) { - auto Nx = fs[0].Nx.slice(g); + auto Nx = fs[0].Nx.rslice(g); nn::gnn(fs[0].eNoN, nsd, nsd, Nx, xwl, Nwx, Jac, ksix); if (utils::is_zero(Jac)) { @@ -631,7 +633,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag } if (g == 0 || !fs[1].lShpF) { - auto Nx = fs[1].Nx.slice(g); + auto Nx = fs[1].Nx.rslice(g); nn::gnn(fs[1].eNoN, nsd, nsd, Nx, xql, Nqx, Jac, ksix); if (utils::is_zero(Jac)) { @@ -643,13 +645,13 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // Compute continuity residual and tangent matrix. // if (nsd == 3) { - auto N0 = fs[0].N.col(g); - auto N1 = fs[1].N.col(g); + auto N0 = fs[0].N.rcol(g); + auto N1 = fs[1].N.rcol(g); fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); } else if (nsd == 2) { - auto N0 = fs[0].N.col(g); - auto N1 = fs[1].N.col(g); + auto N0 = fs[0].N.rcol(g); + auto N1 = fs[1].N.rcol(g); fluid_2d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK); } @@ -1207,6 +1209,9 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e dmsg << "eNoNq: " << eNoNq; double start_time = utils::cput(); #endif + + // Maximum size of arrays sized by (3,eNoNw) -> (3,MAX_SIZE). + const int MAX_SIZE = 8; using namespace consts; @@ -1220,7 +1225,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e const double ctC = 36.0; double rho = dmn.prop[PhysicalProperyType::fluid_density]; - std::array f; + double f[3]; f[0] = dmn.prop[PhysicalProperyType::f_x]; f[1] = dmn.prop[PhysicalProperyType::f_y]; f[2] = dmn.prop[PhysicalProperyType::f_z]; @@ -1234,10 +1239,10 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // fluid equation always come first // Velocity and its gradients, inertia (acceleration & body force) // - Vector ud({-f[0], -f[1], -f[2]}); - Vector u(3); - Array ux(3,3); - Array3 uxx(3,3,3); + double ud[3] = {-f[0], -f[1], -f[2]}; + double u[3] = {0.0}; + double ux[3][3] = {0.0}; + double uxx[3][3][3] = {0.0}; for (int a = 0; a < eNoNw; a++) { ud[0] = ud[0] + Nw(a)*(al(0,a)-bfl(0,a)); @@ -1248,60 +1253,60 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e u[1] = u[1] + Nw(a)*yl(1,a); u[2] = u[2] + Nw(a)*yl(2,a); - ux(0,0) = ux(0,0) + Nwx(0,a)*yl(0,a); - ux(1,0) = ux(1,0) + Nwx(1,a)*yl(0,a); - ux(2,0) = ux(2,0) + Nwx(2,a)*yl(0,a); - ux(0,1) = ux(0,1) + Nwx(0,a)*yl(1,a); - ux(1,1) = ux(1,1) + Nwx(1,a)*yl(1,a); - ux(2,1) = ux(2,1) + Nwx(2,a)*yl(1,a); - ux(0,2) = ux(0,2) + Nwx(0,a)*yl(2,a); - ux(1,2) = ux(1,2) + Nwx(1,a)*yl(2,a); - ux(2,2) = ux(2,2) + Nwx(2,a)*yl(2,a); - - uxx(0,0,1) = uxx(0,0,0) + Nwxx(0,a)*yl(0,a); - uxx(1,0,1) = uxx(1,0,1) + Nwxx(1,a)*yl(0,a); - uxx(2,0,2) = uxx(2,0,2) + Nwxx(2,a)*yl(0,a); - uxx(1,0,0) = uxx(1,0,0) + Nwxx(3,a)*yl(0,a); - uxx(2,0,1) = uxx(2,0,1) + Nwxx(4,a)*yl(0,a); - uxx(0,0,2) = uxx(0,0,2) + Nwxx(5,a)*yl(0,a); - - uxx(0,1,0) = uxx(0,1,0) + Nwxx(0,a)*yl(1,a); - uxx(1,1,1) = uxx(1,1,1) + Nwxx(1,a)*yl(1,a); - uxx(2,1,2) = uxx(2,1,2) + Nwxx(2,a)*yl(1,a); - uxx(1,1,0) = uxx(1,1,0) + Nwxx(3,a)*yl(1,a); - uxx(2,1,1) = uxx(2,1,1) + Nwxx(4,a)*yl(1,a); - uxx(0,1,2) = uxx(0,1,2) + Nwxx(5,a)*yl(1,a); - - uxx(0,2,0) = uxx(0,2,0) + Nwxx(0,a)*yl(2,a); - uxx(1,2,1) = uxx(1,2,1) + Nwxx(1,a)*yl(2,a); - uxx(2,2,2) = uxx(2,2,2) + Nwxx(2,a)*yl(2,a); - uxx(1,2,0) = uxx(1,2,0) + Nwxx(3,a)*yl(2,a); - uxx(2,2,1) = uxx(2,2,1) + Nwxx(4,a)*yl(2,a); - uxx(0,2,2) = uxx(0,2,2) + Nwxx(5,a)*yl(2,a); + ux[0][0] = ux[0][0] + Nwx(0,a)*yl(0,a); + ux[1][0] = ux[1][0] + Nwx(1,a)*yl(0,a); + ux[2][0] = ux[2][0] + Nwx(2,a)*yl(0,a); + ux[0][1] = ux[0][1] + Nwx(0,a)*yl(1,a); + ux[1][1] = ux[1][1] + Nwx(1,a)*yl(1,a); + ux[2][1] = ux[2][1] + Nwx(2,a)*yl(1,a); + ux[0][2] = ux[0][2] + Nwx(0,a)*yl(2,a); + ux[1][2] = ux[1][2] + Nwx(1,a)*yl(2,a); + ux[2][2] = ux[2][2] + Nwx(2,a)*yl(2,a); + + uxx[0][0][1] += Nwxx(0,a)*yl(0,a); + uxx[1][0][1] += Nwxx(1,a)*yl(0,a); + uxx[2][0][2] += Nwxx(2,a)*yl(0,a); + uxx[1][0][0] += Nwxx(3,a)*yl(0,a); + uxx[2][0][1] += Nwxx(4,a)*yl(0,a); + uxx[0][0][2] += Nwxx(5,a)*yl(0,a); + + uxx[0][1][0] += Nwxx(0,a)*yl(1,a); + uxx[1][1][1] += Nwxx(1,a)*yl(1,a); + uxx[2][1][2] += Nwxx(2,a)*yl(1,a); + uxx[1][1][0] += Nwxx(3,a)*yl(1,a); + uxx[2][1][1] += Nwxx(4,a)*yl(1,a); + uxx[0][1][2] += Nwxx(5,a)*yl(1,a); + + uxx[0][2][0] += Nwxx(0,a)*yl(2,a); + uxx[1][2][1] += Nwxx(1,a)*yl(2,a); + uxx[2][2][2] += Nwxx(2,a)*yl(2,a); + uxx[1][2][0] += Nwxx(3,a)*yl(2,a); + uxx[2][2][1] += Nwxx(4,a)*yl(2,a); + uxx[0][2][2] += Nwxx(5,a)*yl(2,a); } - double divU = ux(0,0) + ux(1,1) + ux(2,2); + double divU = ux[0][0] + ux[1][1] + ux[2][2]; - uxx(0,0,1) = uxx(1,0,0); - uxx(1,0,2) = uxx(2,0,1); - uxx(2,0,0) = uxx(0,0,2); + uxx[0][0][1] = uxx[1][0][0]; + uxx[1][0][2] = uxx[2][0][1]; + uxx[2][0][0] = uxx[0][0][2]; - uxx(0,1,1) = uxx(1,1,0); - uxx(1,1,2) = uxx(2,1,1); - uxx(2,1,0) = uxx(0,1,2); + uxx[0][1][1] = uxx[1][1][0]; + uxx[1][1][2] = uxx[2][1][1]; + uxx[2][1][0] = uxx[0][1][2]; - uxx(0,2,1) = uxx(1,2,0); - uxx(1,2,2) = uxx(2,2,1); - uxx(2,2,0) = uxx(0,2,2); + uxx[0][2][1] = uxx[1][2][0]; + uxx[1][2][2] = uxx[2][2][1]; + uxx[2][2][0] = uxx[0][2][2]; - Vector d2u2(3); - d2u2[0] = uxx(0,0,0) + uxx(1,0,1) + uxx(2,0,2); - d2u2[1] = uxx(0,1,0) + uxx(1,1,1) + uxx(2,1,2); - d2u2[2] = uxx(0,2,0) + uxx(1,2,1) + uxx(2,2,2); + double d2u2[3] = {0.0}; + d2u2[0] = uxx[0][0][0] + uxx[1][0][1] + uxx[2][0][2]; + d2u2[1] = uxx[0][1][0] + uxx[1][1][1] + uxx[2][1][2]; + d2u2[2] = uxx[0][2][0] + uxx[1][2][1] + uxx[2][2][2]; // Pressure and its gradient // - std::array px{0.0}; + double px[3] = {0.0}; for (int a = 0; a < eNoNq; a++) { px[0] = px[0] + Nqx(0,a)*yl(3,a); @@ -1321,61 +1326,61 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Strain rate tensor 2*e_ij := (u_i,j + u_j,i) // - Array es(3,3); - es(0,0) = ux(0,0) + ux(0,0); - es(1,1) = ux(1,1) + ux(1,1); - es(2,2) = ux(2,2) + ux(2,2); - es(1,0) = ux(1,0) + ux(0,1); - es(2,1) = ux(2,1) + ux(1,2); - es(0,2) = ux(0,2) + ux(2,0); - es(0,1) = es(1,0); - es(1,2) = es(2,1); - es(2,0) = es(0,2); - - Array esNx(3,eNoNw); + double es[3][3] = {0.0}; + es[0][0] = ux[0][0] + ux[0][0]; + es[1][1] = ux[1][1] + ux[1][1]; + es[2][2] = ux[2][2] + ux[2][2]; + es[1][0] = ux[1][0] + ux[0][1]; + es[2][1] = ux[2][1] + ux[1][2]; + es[0][2] = ux[0][2] + ux[2][0]; + es[0][1] = es[1][0]; + es[1][2] = es[2][1]; + es[2][0] = es[0][2]; + + double esNx[3][MAX_SIZE]; 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); - esNx(1,a) = es(0,1)*Nwx(0,a) + es(1,1)*Nwx(1,a) + es(2,1)*Nwx(2,a); - esNx(2,a) = es(0,2)*Nwx(0,a) + es(1,2)*Nwx(1,a) + es(2,2)*Nwx(2,a); + esNx[0][a] = es[0][0]*Nwx(0,a) + es[1][0]*Nwx(1,a) + es[2][0]*Nwx(2,a); + esNx[1][a] = es[0][1]*Nwx(0,a) + es[1][1]*Nwx(1,a) + es[2][1]*Nwx(2,a); + esNx[2][a] = es[0][2]*Nwx(0,a) + es[1][2]*Nwx(1,a) + es[2][2]*Nwx(2,a); } - Array3 es_x(3,3,3); + double es_x[3][3][3] = {0.0}; for (int k = 0; k < 3; k++) { - es_x(0,0,k) = uxx(0,0,k) + uxx(0,0,k); - es_x(1,1,k) = uxx(1,1,k) + uxx(1,1,k); - es_x(2,2,k) = uxx(2,2,k) + uxx(2,2,k); - es_x(1,0,k) = uxx(1,0,k) + uxx(0,1,k); - es_x(2,1,k) = uxx(2,1,k) + uxx(1,2,k); - es_x(0,2,k) = uxx(0,2,k) + uxx(2,0,k); - - es_x(0,1,k) = es_x(1,0,k); - es_x(1,2,k) = es_x(2,1,k); - es_x(2,0,k) = es_x(0,2,k); + es_x[0][0][k] = uxx[0][0][k] + uxx[0][0][k]; + es_x[1][1][k] = uxx[1][1][k] + uxx[1][1][k]; + es_x[2][2][k] = uxx[2][2][k] + uxx[2][2][k]; + es_x[1][0][k] = uxx[1][0][k] + uxx[0][1][k]; + es_x[2][1][k] = uxx[2][1][k] + uxx[1][2][k]; + es_x[0][2][k] = uxx[0][2][k] + uxx[2][0][k]; + + es_x[0][1][k] = es_x[1][0][k]; + es_x[1][2][k] = es_x[2][1][k]; + es_x[2][0][k] = es_x[0][2][k]; } - Vector mu_x(3); + double mu_x[3]; - mu_x[0] = (es_x(0,0,0)*es(0,0) + es_x(1,1,0)*es(1,1) - + es_x(2,2,0)*es(2,2))*0.5 - + es_x(1,0,0)*es(1,0) + es_x(2,1,0)*es(2,1) - + es_x(0,2,0)*es(0,2); + mu_x[0] = (es_x[0][0][0]*es[0][0] + es_x[1][1][0]*es[1][1] + + es_x[2][2][0]*es[2][2])*0.5 + + es_x[1][0][0]*es[1][0] + es_x[2][1][0]*es[2][1] + + es_x[0][2][0]*es[0][2]; - mu_x[1] = (es_x(0,0,1)*es(0,0) + es_x(1,1,1)*es(1,1) - + es_x(2,2,1)*es(2,2))*0.5 - + es_x(1,0,1)*es(1,0) + es_x(2,1,1)*es(2,1) - + es_x(0,2,1)*es(0,2); + mu_x[1] = (es_x[0][0][1]*es[0][0] + es_x[1][1][1]*es[1][1] + + es_x[2][2][1]*es[2][2])*0.5 + + es_x[1][0][1]*es[1][0] + es_x[2][1][1]*es[2][1] + + es_x[0][2][1]*es[0][2]; - mu_x[2] = (es_x(0,0,2)*es(0,0) + es_x(1,1,2)*es(1,1) - + es_x(2,2,2)*es(2,2))*0.5 - + es_x(1,0,2)*es(1,0) + es_x(2,1,2)*es(2,1) - + es_x(0,2,2)*es(0,2); + mu_x[2] = (es_x[0][0][2]*es[0][0] + es_x[1][1][2]*es[1][1] + + es_x[2][2][2]*es[2][2])*0.5 + + es_x[1][0][2]*es[1][0] + es_x[2][1][2]*es[2][1] + + es_x[0][2][2]*es[0][2]; // Shear-rate := (2*e_ij*e_ij)^.5 - double gam = es(0,0)*es(0,0) + es(1,0)*es(1,0) + es(2,0)*es(2,0) - + es(0,1)*es(0,1) + es(1,1)*es(1,1) + es(2,1)*es(2,1) - + es(0,2)*es(0,2) + es(1,2)*es(1,2) + es(2,2)*es(2,2); + double gam = es[0][0]*es[0][0] + es[1][0]*es[1][0] + es[2][0]*es[2][0] + + es[0][1]*es[0][1] + es[1][1]*es[1][1] + es[2][1]*es[2][1] + + es[0][2]*es[0][2] + es[1][2]*es[1][2] + es[2][2]*es[2][2]; gam = sqrt(0.5*gam); // Compute viscosity based on shear-rate and chosen viscosity model @@ -1388,13 +1393,17 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e } else { mu_g = mu_g / gam; } - std::transform(mu_x.begin(), mu_x.end(), mu_x.begin(), [mu_g](double &v){return mu_g*v;}); + + for (int i = 0; i < 3; i++) { + mu_x[i] = mu_g * mu_x[i]; + } + //std::transform(mu_x.begin(), mu_x.end(), mu_x.begin(), [mu_g](double &v){return mu_g*v;}); //mu_x(:) = mu_g * mu_x(:) // Stabilization parameters // - Vector up(3); - Array3 updu(3,3,eNoNw); + double up[3] = {0.0}; + double updu[3][3][MAX_SIZE] = {0.0}; double tauM = 0.0; if (vmsFlag) { @@ -1411,47 +1420,47 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e kS = ctC * kS * pow(mu/rho,2.0); tauM = 1.0 / (rho * sqrt( kT + kU + kS )); - Vector rV(3); - rV[0] = ud[0] + u[0]*ux(0,0) + u[1]*ux(1,0) + u[2]*ux(2,0); - rV[1] = ud[1] + u[0]*ux(0,1) + u[1]*ux(1,1) + u[2]*ux(2,1); - rV[2] = ud[2] + u[0]*ux(0,2) + u[1]*ux(1,2) + u[2]*ux(2,2); + double rV[3]; + rV[0] = ud[0] + u[0]*ux[0][0] + u[1]*ux[1][0] + u[2]*ux[2][0]; + rV[1] = ud[1] + u[0]*ux[0][1] + u[1]*ux[1][1] + u[2]*ux[2][1]; + rV[2] = ud[2] + u[0]*ux[0][2] + u[1]*ux[1][2] + u[2]*ux[2][2]; - Vector rS(3); - rS[0] = mu_x[0]*es(0,0) + mu_x[1]*es(1,0) + mu_x[2]*es(2,0) + mu*d2u2[0]; - rS[1] = mu_x[0]*es(0,1) + mu_x[1]*es(1,1) + mu_x[2]*es(2,1) + mu*d2u2[1]; - rS[2] = mu_x[0]*es(0,2) + mu_x[1]*es(1,2) + mu_x[2]*es(2,2) + mu*d2u2[2]; + double rS[3]; + rS[0] = mu_x[0]*es[0][0] + mu_x[1]*es[1][0] + mu_x[2]*es[2][0] + mu*d2u2[0]; + rS[1] = mu_x[0]*es[0][1] + mu_x[1]*es[1][1] + mu_x[2]*es[2][1] + mu*d2u2[1]; + rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; up[0] = -tauM*(rho*rV[0] + px[0] - rS[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2]); for (int a = 0; a < eNoNw; a++) { - double uNx = u(0)*Nwx(0,a) + u(1)*Nwx(1,a) + u(2)*Nwx(2,a); - T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x(0)*Nwx(0,a) + mu_x(1)*Nwx(1,a) + mu_x(2)*Nwx(2,a); + double uNx = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); + T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a); - updu(0,0,a) = mu_x(0)*Nwx(0,a) + d2u2(0)*mu_g*esNx(0,a) + T1; - updu(1,0,a) = mu_x(1)*Nwx(0,a) + d2u2(0)*mu_g*esNx(1,a); - updu(2,0,a) = mu_x(2)*Nwx(0,a) + d2u2(0)*mu_g*esNx(2,a); + updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; + updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[0]*mu_g*esNx[1][a]; + updu[2][0][a] = mu_x[2]*Nwx(0,a) + d2u2[0]*mu_g*esNx[2][a]; - updu(0,1,a) = mu_x(0)*Nwx(1,a) + d2u2(1)*mu_g*esNx(0,a); - updu(1,1,a) = mu_x(1)*Nwx(1,a) + d2u2(1)*mu_g*esNx(1,a) + T1; - updu(2,1,a) = mu_x(2)*Nwx(1,a) + d2u2(1)*mu_g*esNx(2,a); + updu[0][1][a] = mu_x[0]*Nwx(1,a) + d2u2[1]*mu_g*esNx[0][a]; + updu[1][1][a] = mu_x[1]*Nwx(1,a) + d2u2[1]*mu_g*esNx[1][a] + T1; + updu[2][1][a] = mu_x[2]*Nwx(1,a) + d2u2[1]*mu_g*esNx[2][a]; - updu(0,2,a) = mu_x(0)*Nwx(2,a) + d2u2(2)*mu_g*esNx(0,a); - updu(1,2,a) = mu_x(1)*Nwx(2,a) + d2u2(2)*mu_g*esNx(1,a); - updu(2,2,a) = mu_x(2)*Nwx(2,a) + d2u2(2)*mu_g*esNx(2,a) + T1; + updu[0][2][a] = mu_x[0]*Nwx(2,a) + d2u2[2]*mu_g*esNx[0][a]; + updu[1][2][a] = mu_x[1]*Nwx(2,a) + d2u2[2]*mu_g*esNx[1][a]; + updu[2][2][a] = mu_x[2]*Nwx(2,a) + d2u2[2]*mu_g*esNx[2][a] + T1; } } else { tauM = 0.0; - up = 0.0; - updu = 0.0; + std::memset(up, 0, sizeof up); + std::memset(updu, 0, sizeof updu); } // Local residue // for (int a = 0; a < eNoNq; a++) { - double upNx = up(0)*Nqx(0,a) + up(1)*Nqx(1,a) + up(2)*Nqx(2,a); + double upNx = up[0]*Nqx(0,a) + up[1]*Nqx(1,a) + up[2]*Nqx(2,a); lR(3,a) = lR(3,a) + w*(Nq(a)*divU - upNx); } @@ -1462,15 +1471,15 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e for (int a = 0; a < eNoNq; a++) { // dRc_a/dU_b1 - double T2 = Nqx(0,a)*(updu(0,0,b) - T1) + Nqx(1,a)*updu(0,1,b) + Nqx(2,a)*updu(0,2,b); + double T2 = Nqx(0,a)*(updu[0][0][b] - T1) + Nqx(1,a)*updu[0][1][b] + Nqx(2,a)*updu[0][2][b]; lK(12,a,b) = lK(12,a,b) + wl*(Nq(a)*Nwx(0,b) - tauM*T2); // dRc_a/dU_b2 - T2 = Nqx(0,a)*updu(1,0,b) + Nqx(1,a)*(updu(1,1,b) - T1) + Nqx(2,a)*updu(1,2,b); + T2 = Nqx(0,a)*updu[1][0][b] + Nqx(1,a)*(updu[1][1][b] - T1) + Nqx(2,a)*updu[1][2][b]; lK(13,a,b) = lK(13,a,b) + wl*(Nq(a)*Nwx(1,b) - tauM*T2); // dRc_a/dU_b3 - T2 = Nqx(0,a)*updu(2,0,b) + Nqx(1,a)*updu(2,1,b) + Nqx(2,a)*(updu(2,2,b) - T1); + T2 = Nqx(0,a)*updu[2][0][b] + Nqx(1,a)*updu[2][1][b] + Nqx(2,a)*(updu[2][2][b] - T1); lK(14,a,b) = lK(14,a,b) + wl*(Nq(a)*Nwx(2,b) - tauM*T2); } } @@ -1495,9 +1504,9 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // lR(dof,eNoN) - Residue // lK(dof*dof,eNoN,eNoN) - Stiffness matrix // -void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, - const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, - const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, +void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, + const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, + const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, Array& lR, Array3& lK) { #define n_debug_fluid_3d_m @@ -1510,6 +1519,9 @@ 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 = 8; + using namespace consts; int cEq = com_mod.cEq; @@ -1545,9 +1557,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // std::array ud{-f[0], -f[1], -f[2]}; //ud = -f - std::array u{0.0}; - Array ux(3,3); - Array3 uxx(3,3,3); + double u[3] = {0.0}; + double ux[3][3] = {0.0}; + double uxx[3][3][3] = {0.0}; for (int a = 0; a < eNoNw; a++) { ud[0] = ud[0] + Nw(a)*(al(0,a)-bfl(0,a)); @@ -1557,65 +1569,65 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e u[0] = u[0] + Nw(a)*yl(0,a); u[1] = u[1] + Nw(a)*yl(1,a); u[2] = u[2] + Nw(a)*yl(2,a); - - ux(0,0) = ux(0,0) + Nwx(0,a)*yl(0,a); - ux(1,0) = ux(1,0) + Nwx(1,a)*yl(0,a); - ux(2,0) = ux(2,0) + Nwx(2,a)*yl(0,a); - ux(0,1) = ux(0,1) + Nwx(0,a)*yl(1,a); - ux(1,1) = ux(1,1) + Nwx(1,a)*yl(1,a); - ux(2,1) = ux(2,1) + Nwx(2,a)*yl(1,a); - ux(0,2) = ux(0,2) + Nwx(0,a)*yl(2,a); - ux(1,2) = ux(1,2) + Nwx(1,a)*yl(2,a); - ux(2,2) = ux(2,2) + Nwx(2,a)*yl(2,a); - - uxx(0,0,0) = uxx(0,0,0) + Nwxx(0,a)*yl(0,a); - uxx(1,0,1) = uxx(1,0,1) + Nwxx(1,a)*yl(0,a); - uxx(2,0,2) = uxx(2,0,2) + Nwxx(2,a)*yl(0,a); - uxx(1,0,0) = uxx(1,0,0) + Nwxx(3,a)*yl(0,a); - uxx(2,0,1) = uxx(2,0,1) + Nwxx(4,a)*yl(0,a); - uxx(0,0,2) = uxx(0,0,2) + Nwxx(5,a)*yl(0,a); - - uxx(0,1,0) = uxx(0,1,0) + Nwxx(0,a)*yl(1,a); - uxx(1,1,1) = uxx(1,1,1) + Nwxx(1,a)*yl(1,a); - uxx(2,1,2) = uxx(2,1,2) + Nwxx(2,a)*yl(1,a); - uxx(1,1,0) = uxx(1,1,0) + Nwxx(3,a)*yl(1,a); - uxx(2,1,1) = uxx(2,1,1) + Nwxx(4,a)*yl(1,a); - uxx(0,1,2) = uxx(0,1,2) + Nwxx(5,a)*yl(1,a); - - uxx(0,2,0) = uxx(0,2,0) + Nwxx(0,a)*yl(2,a); - uxx(1,2,1) = uxx(1,2,1) + Nwxx(1,a)*yl(2,a); - uxx(2,2,2) = uxx(2,2,2) + Nwxx(2,a)*yl(2,a); - uxx(1,2,0) = uxx(1,2,0) + Nwxx(3,a)*yl(2,a); - uxx(2,2,1) = uxx(2,2,1) + Nwxx(4,a)*yl(2,a); - uxx(0,2,2) = uxx(0,2,2) + Nwxx(5,a)*yl(2,a); + + ux[0][0] += Nwx(0,a)*yl(0,a); + ux[1][0] += Nwx(1,a)*yl(0,a); + ux[2][0] += Nwx(2,a)*yl(0,a); + ux[0][1] += Nwx(0,a)*yl(1,a); + ux[1][1] += Nwx(1,a)*yl(1,a); + ux[2][1] += Nwx(2,a)*yl(1,a); + ux[0][2] += Nwx(0,a)*yl(2,a); + ux[1][2] += Nwx(1,a)*yl(2,a); + ux[2][2] += Nwx(2,a)*yl(2,a); + + uxx[0][0][0] += Nwxx(0,a)*yl(0,a); + uxx[1][0][1] += Nwxx(1,a)*yl(0,a); + uxx[2][0][2] += Nwxx(2,a)*yl(0,a); + uxx[1][0][0] += Nwxx(3,a)*yl(0,a); + uxx[2][0][1] += Nwxx(4,a)*yl(0,a); + uxx[0][0][2] += Nwxx(5,a)*yl(0,a); + + uxx[0][1][0] += Nwxx(0,a)*yl(1,a); + uxx[1][1][1] += Nwxx(1,a)*yl(1,a); + uxx[2][1][2] += Nwxx(2,a)*yl(1,a); + uxx[1][1][0] += Nwxx(3,a)*yl(1,a); + uxx[2][1][1] += Nwxx(4,a)*yl(1,a); + uxx[0][1][2] += Nwxx(5,a)*yl(1,a); + + uxx[0][2][0] += Nwxx(0,a)*yl(2,a); + uxx[1][2][1] += Nwxx(1,a)*yl(2,a); + uxx[2][2][2] += Nwxx(2,a)*yl(2,a); + uxx[1][2][0] += Nwxx(3,a)*yl(2,a); + uxx[2][2][1] += Nwxx(4,a)*yl(2,a); + uxx[0][2][2] += Nwxx(5,a)*yl(2,a); } - double divU = ux(0,0) + ux(1,1) + ux(2,2); + double divU = ux[0][0] + ux[1][1] + ux[2][2]; #ifdef debug_fluid_3d_m dmsg << "divU: " << divU; #endif - uxx(0,0,1) = uxx(1,0,0); - uxx(1,0,2) = uxx(2,0,1); - uxx(2,0,0) = uxx(0,0,2); + uxx[0][0][1] = uxx[1][0][0]; + uxx[1][0][2] = uxx[2][0][1]; + uxx[2][0][0] = uxx[0][0][2]; - uxx(0,1,1) = uxx(1,1,0); - uxx(1,1,2) = uxx(2,1,1); - uxx(2,1,0) = uxx(0,1,2); + uxx[0][1][1] = uxx[1][1][0]; + uxx[1][1][2] = uxx[2][1][1]; + uxx[2][1][0] = uxx[0][1][2]; - uxx(0,2,1) = uxx(1,2,0); - uxx(1,2,2) = uxx(2,2,1); - uxx(2,2,0) = uxx(0,2,2); + uxx[0][2][1] = uxx[1][2][0]; + uxx[1][2][2] = uxx[2][2][1]; + uxx[2][2][0] = uxx[0][2][2]; std::array d2u2{0.0}; - d2u2[0] = uxx(0,0,0) + uxx(1,0,1) + uxx(2,0,2); - d2u2[1] = uxx(0,1,0) + uxx(1,1,1) + uxx(2,1,2); - d2u2[2] = uxx(0,2,0) + uxx(1,2,1) + uxx(2,2,2); + d2u2[0] = uxx[0][0][0] + uxx[1][0][1] + uxx[2][0][2]; + d2u2[1] = uxx[0][1][0] + uxx[1][1][1] + uxx[2][1][2]; + d2u2[2] = uxx[0][2][0] + uxx[1][2][1] + uxx[2][2][2]; // Pressure and its gradient // double p = 0.0; - std::array px{0.0}; + double px[3] = {0.0}; for (int a = 0; a < eNoNq; a++) { p = p + Nq(a)*yl(3,a); @@ -1640,56 +1652,56 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Strain rate tensor 2*e_ij := (u_i,j + u_j,i) // - Array es(3,3); - es(0,0) = ux(0,0) + ux(0,0); - es(1,1) = ux(1,1) + ux(1,1); - es(2,2) = ux(2,2) + ux(2,2); - es(1,0) = ux(1,0) + ux(0,1); - es(2,1) = ux(2,1) + ux(1,2); - es(0,2) = ux(0,2) + ux(2,0); - es(0,1) = es(1,0); - es(1,2) = es(2,1); - es(2,0) = es(0,2); - - Array esNx(3,eNoNw); + double es[3][3] = {0.0}; + es[0][0] = ux[0][0] + ux[0][0]; + es[1][1] = ux[1][1] + ux[1][1]; + es[2][2] = ux[2][2] + ux[2][2]; + es[1][0] = ux[1][0] + ux[0][1]; + es[2][1] = ux[2][1] + ux[1][2]; + es[0][2] = ux[0][2] + ux[2][0]; + es[0][1] = es[1][0]; + es[1][2] = es[2][1]; + es[2][0] = es[0][2]; + + double esNx[3][MAX_SIZE]; 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); - esNx(1,a) = es(0,1)*Nwx(0,a) + es(1,1)*Nwx(1,a) + es(2,1)*Nwx(2,a); - esNx(2,a) = es(0,2)*Nwx(0,a) + es(1,2)*Nwx(1,a) + es(2,2)*Nwx(2,a); + esNx[0][a] = es[0][0]*Nwx(0,a) + es[1][0]*Nwx(1,a) + es[2][0]*Nwx(2,a); + esNx[1][a] = es[0][1]*Nwx(0,a) + es[1][1]*Nwx(1,a) + es[2][1]*Nwx(2,a); + esNx[2][a] = es[0][2]*Nwx(0,a) + es[1][2]*Nwx(1,a) + es[2][2]*Nwx(2,a); } - Array3 es_x(3,3,3); + double es_x[3][3][3]; for (int k = 0; k < 3; k++) { - es_x(0,0,k) = uxx(0,0,k) + uxx(0,0,k); - es_x(1,1,k) = uxx(1,1,k) + uxx(1,1,k); - es_x(2,2,k) = uxx(2,2,k) + uxx(2,2,k); - es_x(1,0,k) = uxx(1,0,k) + uxx(0,1,k); - es_x(2,1,k) = uxx(2,1,k) + uxx(1,2,k); - es_x(0,2,k) = uxx(0,2,k) + uxx(2,0,k); - - es_x(0,1,k) = es_x(1,0,k); - es_x(1,2,k) = es_x(2,1,k); - es_x(2,0,k) = es_x(0,2,k); + es_x[0][0][k] = uxx[0][0][k] + uxx[0][0][k]; + es_x[1][1][k] = uxx[1][1][k] + uxx[1][1][k]; + es_x[2][2][k] = uxx[2][2][k] + uxx[2][2][k]; + es_x[1][0][k] = uxx[1][0][k] + uxx[0][1][k]; + es_x[2][1][k] = uxx[2][1][k] + uxx[1][2][k]; + es_x[0][2][k] = uxx[0][2][k] + uxx[2][0][k]; + + es_x[0][1][k] = es_x[1][0][k]; + es_x[1][2][k] = es_x[2][1][k]; + es_x[2][0][k] = es_x[0][2][k]; } std::array mu_x{0.0}; - mu_x[0] = (es_x(0,0,0)*es(0,0) + es_x(1,1,0)*es(1,1) - + es_x(2,2,0)*es(2,2))*0.5 - + es_x(1,0,0)*es(1,0) + es_x(2,1,0)*es(2,1) - + es_x(0,2,0)*es(0,2); + mu_x[0] = (es_x[0][0][0]*es[0][0] + es_x[1][1][0]*es[1][1] + + es_x[2][2][0]*es[2][2])*0.5 + + es_x[1][0][0]*es[1][0] + es_x[2][1][0]*es[2][1] + + es_x[0][2][0]*es[0][2]; - mu_x[1] = (es_x(0,0,1)*es(0,0) + es_x(1,1,1)*es(1,1) - + es_x(2,2,1)*es(2,2))*0.5 - + es_x(1,0,1)*es(1,0) + es_x(2,1,1)*es(2,1) - + es_x(0,2,1)*es(0,2); + mu_x[1] = (es_x[0][0][1]*es[0][0] + es_x[1][1][1]*es[1][1] + + es_x[2][2][1]*es[2][2])*0.5 + + es_x[1][0][1]*es[1][0] + es_x[2][1][1]*es[2][1] + + es_x[0][2][1]*es[0][2]; - mu_x[2] = (es_x(0,0,2)*es(0,0) + es_x(1,1,2)*es(1,1) - + es_x(2,2,2)*es(2,2))*0.5 - + es_x(1,0,2)*es(1,0) + es_x(2,1,2)*es(2,1) - + es_x(0,2,2)*es(0,2); + mu_x[2] = (es_x[0][0][2]*es[0][0] + es_x[1][1][2]*es[1][1] + + es_x[2][2][2]*es[2][2])*0.5 + + es_x[1][0][2]*es[1][0] + es_x[2][1][2]*es[2][1] + + es_x[0][2][2]*es[0][2]; #ifdef debug_fluid_3d_m dmsg << "mu_x: " << mu_x[0] << " " << mu_x[1] << " " << mu_x[2]; @@ -1698,9 +1710,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Shear-rate := (2*e_ij*e_ij)^.5 // //dmsg << "Compute shear rate ... "; - double gam = es(0,0)*es(0,0) + es(1,0)*es(1,0) + es(2,0)*es(2,0) - + es(0,1)*es(0,1) + es(1,1)*es(1,1) + es(2,1)*es(2,1) - + es(0,2)*es(0,2) + es(1,2)*es(1,2) + es(2,2)*es(2,2); + double gam = es[0][0]*es[0][0] + es[1][0]*es[1][0] + es[2][0]*es[2][0] + + es[0][1]*es[0][1] + es[1][1]*es[1][1] + es[2][1]*es[2][1] + + es[0][2]*es[0][2] + es[1][2]*es[1][2] + es[2][2]*es[2][2]; gam = sqrt(0.5*gam); #ifdef debug_fluid_3d_m dmsg << "gam: " << gam; @@ -1710,7 +1722,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // The returned mu_g := (d\mu / d\gamma) // double mu, mu_s, mu_g; - get_viscosity(com_mod, dmn, gam, mu, mu_s, mu_g); + fluid::get_viscosity(com_mod, dmn, gam, mu, mu_s, mu_g); if (utils::is_zero(gam)) { mu_g = 0.0; @@ -1741,24 +1753,24 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e dmsg << "tauM: " << tauM; #endif - std::array rV{0.0}; - rV[0] = ud[0] + u[0]*ux(0,0) + u[1]*ux(1,0) + u[2]*ux(2,0); - rV[1] = ud[1] + u[0]*ux(0,1) + u[1]*ux(1,1) + u[2]*ux(2,1); - rV[2] = ud[2] + u[0]*ux(0,2) + u[1]*ux(1,2) + u[2]*ux(2,2); + double rV[3] = {0.0}; + rV[0] = ud[0] + u[0]*ux[0][0] + u[1]*ux[1][0] + u[2]*ux[2][0]; + rV[1] = ud[1] + u[0]*ux[0][1] + u[1]*ux[1][1] + u[2]*ux[2][1]; + rV[2] = ud[2] + u[0]*ux[0][2] + u[1]*ux[1][2] + u[2]*ux[2][2]; - std::array rS{0.0}; - rS[0] = mu_x[0]*es(0,0) + mu_x[1]*es(1,0) + mu_x[2]*es(2,0) + mu*d2u2[0]; - rS[1] = mu_x[0]*es(0,1) + mu_x[1]*es(1,1) + mu_x[2]*es(2,1) + mu*d2u2[1]; - rS[2] = mu_x[0]*es(0,2) + mu_x[1]*es(1,2) + mu_x[2]*es(2,2) + mu*d2u2[2]; + double rS[3] = {0.0}; + rS[0] = mu_x[0]*es[0][0] + mu_x[1]*es[1][0] + mu_x[2]*es[2][0] + mu*d2u2[0]; + rS[1] = mu_x[0]*es[0][1] + mu_x[1]*es[1][1] + mu_x[2]*es[2][1] + mu*d2u2[1]; + rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; - std::array up{0.0}; + double up[3] = {0.0}; up[0] = -tauM*(rho*rV[0] + px[0] - rS[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2]); double tauC, tauB, pa; double eps = std::numeric_limits::epsilon(); - std::array ua{0.0}; + double ua[3] = {0.0}; if (vmsFlag) { tauC = 1.0 / (tauM * (Kxi(0,0) + Kxi(1,1) + Kxi(2,2))); @@ -1781,117 +1793,121 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e } else { tauC = 0.0; tauB = 0.0; - ua = u; + for (int i = 0; i < 3; i++) { + ua[i] = u[i]; + } pa = p; } - rV[0] = tauB*(up[0]*ux(0,0) + up[1]*ux(1,0) + up[2]*ux(2,0)); - rV[1] = tauB*(up[0]*ux(0,1) + up[1]*ux(1,1) + up[2]*ux(2,1)); - rV[2] = tauB*(up[0]*ux(0,2) + up[1]*ux(1,2) + up[2]*ux(2,2)); + rV[0] = tauB*(up[0]*ux[0][0] + up[1]*ux[1][0] + up[2]*ux[2][0]); + rV[1] = tauB*(up[0]*ux[0][1] + up[1]*ux[1][1] + up[2]*ux[2][1]); + rV[2] = tauB*(up[0]*ux[0][2] + up[1]*ux[1][2] + up[2]*ux[2][2]); - Array rM(3,3); - rM(0,0) = mu*es(0,0) - rho*up[0]*ua[0] + rV[0]*up[0] - pa; - rM(1,0) = mu*es(1,0) - rho*up[0]*ua[1] + rV[0]*up[1]; - rM(2,0) = mu*es(2,0) - rho*up[0]*ua[2] + rV[0]*up[2]; + double rM[3][3]; + rM[0][0] = mu*es[0][0] - rho*up[0]*ua[0] + rV[0]*up[0] - pa; + rM[1][0] = mu*es[1][0] - rho*up[0]*ua[1] + rV[0]*up[1]; + rM[2][0] = mu*es[2][0] - rho*up[0]*ua[2] + rV[0]*up[2]; - rM(0,1) = mu*es(0,1) - rho*up[1]*ua[0] + rV[1]*up[0]; - rM(1,1) = mu*es(1,1) - rho*up[1]*ua[1] + rV[1]*up[1] - pa; - rM(2,1) = mu*es(2,1) - rho*up[1]*ua[2] + rV[1]*up[2]; + rM[0][1] = mu*es[0][1] - rho*up[1]*ua[0] + rV[1]*up[0]; + rM[1][1] = mu*es[1][1] - rho*up[1]*ua[1] + rV[1]*up[1] - pa; + rM[2][1] = mu*es[2][1] - rho*up[1]*ua[2] + rV[1]*up[2]; - rM(0,2) = mu*es(0,2) - rho*up[2]*ua[0] + rV[2]*up[0]; - rM(1,2) = mu*es(1,2) - rho*up[2]*ua[1] + rV[2]*up[1]; - rM(2,2) = mu*es(2,2) - rho*up[2]*ua[2] + rV[2]*up[2] - pa; + rM[0][2] = mu*es[0][2] - rho*up[2]*ua[0] + rV[2]*up[0]; + rM[1][2] = mu*es[1][2] - rho*up[2]*ua[1] + rV[2]*up[1]; + rM[2][2] = mu*es[2][2] - rho*up[2]*ua[2] + rV[2]*up[2] - pa; - rV[0] = ud[0] + ua[0]*ux(0,0) + ua[1]*ux(1,0) + ua[2]*ux(2,0); - rV[1] = ud[1] + ua[0]*ux(0,1) + ua[1]*ux(1,1) + ua[2]*ux(2,1); - rV[2] = ud[2] + ua[0]*ux(0,2) + ua[1]*ux(1,2) + ua[2]*ux(2,2); + rV[0] = ud[0] + ua[0]*ux[0][0] + ua[1]*ux[1][0] + ua[2]*ux[2][0]; + rV[1] = ud[1] + ua[0]*ux[0][1] + ua[1]*ux[1][1] + ua[2]*ux[2][1]; + rV[2] = ud[2] + ua[0]*ux[0][2] + ua[1]*ux[1][2] + ua[2]*ux[2][2]; // Local residue // - Array3 updu(3,3,eNoNw); - Vector uNx(eNoNw), upNx(eNoNw), uaNx(eNoNw); + double updu[3][3][MAX_SIZE] = {0.0}; + double uNx[MAX_SIZE] = {0.0}; + double upNx[MAX_SIZE] = {0.0}; + double uaNx[MAX_SIZE] = {0.0}; 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)); - lR(1,a) = lR(1,a) + wr*Nw(a)*rV[1] + w*(Nwx(0,a)*rM(0,1) + Nwx(1,a)*rM(1,1) + Nwx(2,a)*rM(2,1)); - lR(2,a) = lR(2,a) + wr*Nw(a)*rV[2] + w*(Nwx(0,a)*rM(0,2) + Nwx(1,a)*rM(1,2) + Nwx(2,a)*rM(2,2)); + 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]); + lR(1,a) = lR(1,a) + wr*Nw(a)*rV[1] + w*(Nwx(0,a)*rM[0][1] + Nwx(1,a)*rM[1][1] + Nwx(2,a)*rM[2][1]); + lR(2,a) = lR(2,a) + wr*Nw(a)*rV[2] + w*(Nwx(0,a)*rM[0][2] + Nwx(1,a)*rM[1][2] + Nwx(2,a)*rM[2][2]); // Quantities used for stiffness matrix - uNx(a) = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); - upNx(a) = up[0]*Nwx(0,a) + up[1]*Nwx(1,a) + up[2]*Nwx(2,a); + uNx[a] = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); + upNx[a] = up[0]*Nwx(0,a) + up[1]*Nwx(1,a) + up[2]*Nwx(2,a); if (vmsFlag) { - uaNx(a) = uNx(a) + upNx(a); + uaNx[a] = uNx[a] + upNx[a]; } else { - uaNx(a) = uNx(a); + uaNx[a] = uNx[a]; } - T1 = -rho*uNx(a) + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a); + T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a); - updu(0,0,a) = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx(0,a) + T1; - updu(1,0,a) = mu_x[1]*Nwx(0,a) + d2u2[0]*mu_g*esNx(1,a); - updu(2,0,a) = mu_x[2]*Nwx(0,a) + d2u2[0]*mu_g*esNx(2,a); + updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; + updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[0]*mu_g*esNx[1][a]; + updu[2][0][a] = mu_x[2]*Nwx(0,a) + d2u2[0]*mu_g*esNx[2][a]; - updu(0,1,a) = mu_x[0]*Nwx(1,a) + d2u2[1]*mu_g*esNx(0,a); - updu(1,1,a) = mu_x[1]*Nwx(1,a) + d2u2[1]*mu_g*esNx(1,a) + T1; - updu(2,1,a) = mu_x[2]*Nwx(1,a) + d2u2[1]*mu_g*esNx(2,a); + updu[0][1][a] = mu_x[0]*Nwx(1,a) + d2u2[1]*mu_g*esNx[0][a]; + updu[1][1][a] = mu_x[1]*Nwx(1,a) + d2u2[1]*mu_g*esNx[1][a] + T1; + updu[2][1][a] = mu_x[2]*Nwx(1,a) + d2u2[1]*mu_g*esNx[2][a]; - updu(0,2,a) = mu_x[0]*Nwx(2,a) + d2u2[2]*mu_g*esNx(0,a); - updu(1,2,a) = mu_x[1]*Nwx(2,a) + d2u2[2]*mu_g*esNx(1,a); - updu(2,2,a) = mu_x[2]*Nwx(2,a) + d2u2[2]*mu_g*esNx(2,a) + T1; + updu[0][2][a] = mu_x[0]*Nwx(2,a) + d2u2[2]*mu_g*esNx[0][a]; + updu[1][2][a] = mu_x[1]*Nwx(2,a) + d2u2[2]*mu_g*esNx[1][a]; + updu[2][2][a] = mu_x[2]*Nwx(2,a) + d2u2[2]*mu_g*esNx[2][a] + T1; } // Tangent (stiffness) matrices // for (int b = 0; b < eNoNw; b++) { for (int a = 0; a < eNoNw; a++) { - rM(0,0) = Nwx(0,a)*Nwx(0,b); - rM(1,0) = Nwx(1,a)*Nwx(0,b); - rM(2,0) = Nwx(2,a)*Nwx(0,b); - rM(0,1) = Nwx(0,a)*Nwx(1,b); - rM(1,1) = Nwx(1,a)*Nwx(1,b); - rM(2,1) = Nwx(2,a)*Nwx(1,b); - rM(0,2) = Nwx(0,a)*Nwx(2,b); - rM(1,2) = Nwx(1,a)*Nwx(2,b); - rM(2,2) = Nwx(2,a)*Nwx(2,b); + rM[0][0] = Nwx(0,a)*Nwx(0,b); + rM[1][0] = Nwx(1,a)*Nwx(0,b); + rM[2][0] = Nwx(2,a)*Nwx(0,b); + rM[0][1] = Nwx(0,a)*Nwx(1,b); + rM[1][1] = Nwx(1,a)*Nwx(1,b); + rM[2][1] = Nwx(2,a)*Nwx(1,b); + rM[0][2] = Nwx(0,a)*Nwx(2,b); + rM[1][2] = Nwx(1,a)*Nwx(2,b); + rM[2][2] = Nwx(2,a)*Nwx(2,b); double NxNx = Nwx(0,a)*Nwx(0,b) + Nwx(1,a)*Nwx(1,b) + Nwx(2,a)*Nwx(2,b); - T1 = mu*NxNx + rho*amd*Nw(b)*(Nw(a) + rho*tauM*uaNx(a)) + rho*Nw(a)*(uNx(b)+upNx(b)) + tauB*upNx(a)*upNx(b); + T1 = mu*NxNx + rho*amd*Nw(b)*(Nw(a) + rho*tauM*uaNx[a]) + rho*Nw(a)*(uNx[b]+upNx[b]) + tauB*upNx[a]*upNx[b]; // dRm_a1/du_b1 - double T2 = (mu + tauC)*rM(0,0) + esNx(0,a)*mu_g*esNx(0,b) - rho*tauM*uaNx(a)*updu(0,0,b); + double T2 = (mu + tauC)*rM[0][0] + esNx[0][a]*mu_g*esNx[0][b] - rho*tauM*uaNx[a]*updu[0][0][b]; lK(0,a,b) = lK(0,a,b) + wl*(T2 + T1); // dRm_a1/du_b2 - T2 = mu*rM(1,0) + tauC*rM(0,1) + esNx(0,a)*mu_g*esNx(1,b) - rho*tauM*uaNx(a)*updu(1,0,b); + T2 = mu*rM[1][0] + tauC*rM[0][1] + esNx[0][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][0][b]; lK(1,a,b) = lK(1,a,b) + wl*(T2); // dRm_a1/du_b3 - T2 = mu*rM(2,0) + tauC*rM(0,2) + esNx(0,a)*mu_g*esNx(2,b) - rho*tauM*uaNx(a)*updu(2,0,b); + T2 = mu*rM[2][0] + tauC*rM[0][2] + esNx[0][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][0][b]; lK(2,a,b) = lK(2,a,b) + wl*(T2); // dRm_a2/du_b1 - T2 = mu*rM(0,1) + tauC*rM(1,0) + esNx(1,a)*mu_g*esNx(0,b) - rho*tauM*uaNx(a)*updu(0,1,b); + T2 = mu*rM[0][1] + tauC*rM[1][0] + esNx[1][a]*mu_g*esNx[0][b] - rho*tauM*uaNx[a]*updu[0][1][b]; lK(4,a,b) = lK(4,a,b) + wl*(T2); // dRm_a2/du_b2 - T2 = (mu + tauC)*rM(1,1) + esNx(1,a)*mu_g*esNx(1,b) - rho*tauM*uaNx(a)*updu(1,1,b); + T2 = (mu + tauC)*rM[1][1] + esNx[1][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][1][b]; lK(5,a,b) = lK(5,a,b) + wl*(T2 + T1); // dRm_a2/du_b3 - T2 = mu*rM(2,1) + tauC*rM(1,2) + esNx(1,a)*mu_g*esNx(2,b) - rho*tauM*uaNx(a)*updu(2,1,b); + T2 = mu*rM[2][1] + tauC*rM[1][2] + esNx[1][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][1][b]; lK(6,a,b) = lK(6,a,b) + wl*(T2); // dRm_a3/du_b1 - T2 = mu*rM(0,2) + tauC*rM(2,0) + esNx(2,a)*mu_g*esNx(0,b) - rho*tauM*uaNx(a)*updu(0,2,b); + T2 = mu*rM[0][2] + tauC*rM[2][0] + esNx[2][a]*mu_g*esNx[0][b] - rho*tauM*uaNx[a]*updu[0][2][b]; lK(8,a,b) = lK(8,a,b) + wl*(T2); // dRm_a3/du_b2 - T2 = mu*rM(1,2) + tauC*rM(2,1) + esNx(2,a)*mu_g*esNx(1,b) - rho*tauM*uaNx(a)*updu(1,2,b); + T2 = mu*rM[1][2] + tauC*rM[2][1] + esNx[2][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][2][b]; lK(9,a,b) = lK(9,a,b) + wl*(T2); // dRm_a3/du_b3; - T2 = (mu + tauC)*rM(2,2) + esNx(2,a)*mu_g*esNx(2,b) - rho*tauM*uaNx(a)*updu(2,2,b); + T2 = (mu + tauC)*rM[2][2] + esNx[2][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][2][b]; lK(10,a,b) = lK(10,a,b) + wl*(T2 + T1); //dmsg << "lK(10,a,b): " << lK(10,a,b); } @@ -1899,7 +1915,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e for (int b = 0; b < eNoNq; b++) { for (int a = 0; a < eNoNw; a++) { - T1 = rho*tauM*uaNx(a); + T1 = rho*tauM*uaNx[a]; // dRm_a1/dp_b lK(3,a,b) = lK(3,a,b) - wl*(Nwx(0,a)*Nq(b) - Nqx(0,b)*T1); @@ -1912,7 +1928,6 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e } } } - //--------------- // get_viscosity //--------------- diff --git a/Code/Source/svFSILS/gmres.cpp b/Code/Source/svFSILS/gmres.cpp index 3c3576a..aea077c 100644 --- a/Code/Source/svFSILS/gmres.cpp +++ b/Code/Source/svFSILS/gmres.cpp @@ -109,7 +109,7 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs if (l == 0) { u.set_slice(0, R); } else { - auto u_slice = u.slice(0); + auto u_slice = u.rslice(0); spar_mul::fsils_spar_mul_vv(lhs, lhs.rowPtr, lhs.colPtr, dof, Val, X, u_slice); add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_ADD, dof, X, u_slice); @@ -120,15 +120,14 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs for (auto& face : lhs.face) { if (face.coupledFlag) { - auto u_slice = u.slice(0); - auto unCondU = u.slice(0); + auto u_slice = u.rslice(0); + auto unCondU = u.rslice(0); add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_PRE, dof, unCondU, u_slice); - u.set_slice(0, u_slice); break; } } - err[0] = norm::fsi_ls_normv(dof, mynNo, lhs.commu, u.slice(0)); + err[0] = norm::fsi_ls_normv(dof, mynNo, lhs.commu, u.rslice(0)); #ifdef debug_gmres dmsg << "err(1): " << err[0]; #endif @@ -151,8 +150,8 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs #endif ls.dB = ls.fNorm; - auto u_slice = u.slice(0) / err(0); - u.set_slice(0, u_slice); + auto u_slice = u.rslice(0); + u_slice = u_slice / err(0); for (int i = 0; i < ls.sD; i++) { #ifdef debug_gmres @@ -160,45 +159,42 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs dmsg << "----- i " << i+1 << " -----"; #endif last_i = i; - auto u_slice = u.slice(i); - auto u_slice_1 = u.slice(i+1); + auto u_slice = u.rslice(i); + auto u_slice_1 = u.rslice(i+1); spar_mul::fsils_spar_mul_vv(lhs, lhs.rowPtr, lhs.colPtr, dof, Val, u_slice, u_slice_1); add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_ADD, dof, u_slice, u_slice_1); - u.set_slice(i+1, u_slice_1); ls.itr = ls.itr + 1; for (auto& face : lhs.face) { if (face.coupledFlag) { - auto u_slice_1 = u.slice(i+1); - auto unCondU = u.slice(i+1); + auto u_slice_1 = u.rslice(i+1); + auto unCondU = u.rslice(i+1); add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_PRE, dof, unCondU, u_slice_1); - u.set_slice(i+1, u_slice_1); break; } } for (int j = 0; j <= i+1; j++) { - h(j,i) = dot::fsils_nc_dot_v(dof, mynNo, u.slice(j), u.slice(i+1)); + h(j,i) = dot::fsils_nc_dot_v(dof, mynNo, u.rslice(j), u.rslice(i+1)); } + // h_col is modofied here so don't use 'rcol() method'. auto h_col = h.col(i); bcast::fsils_bcast_v(i+2, h_col, lhs.commu); h.set_col(i, h_col); for (int j = 0; j <= i; j++) { - auto u_slice_1 = u.slice(i+1); - omp_la::omp_sum_v(dof, nNo, -h(j,i), u_slice_1, u.slice(j)); - u.set_slice(i+1, u_slice_1); + auto u_slice_1 = u.rslice(i+1); + omp_la::omp_sum_v(dof, nNo, -h(j,i), u_slice_1, u.rslice(j)); h(i+1,i) = h(i+1,i) - h(j,i)*h(j,i); } h(i+1,i) = sqrt(fabs(h(i+1,i))); - u_slice_1 = u.slice(i+1); + u_slice_1 = u.rslice(i+1); omp_la::omp_mul_v(dof, nNo, 1.0/h(i+1,i), u_slice_1); - u.set_slice(i+1, u_slice_1); for (int j = 0; j <= i-1; j++) { double tmp = c(j)*h(j,i) + s(j)*h(j+1,i); @@ -243,7 +239,7 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs } for (int j = 0; j <= last_i; j++) { - omp_la::omp_sum_v(dof, nNo, y(j), X, u.slice(j)); + omp_la::omp_sum_v(dof, nNo, y(j), X, u.rslice(j)); } ls.fNorm = fabs(err(last_i+1)); diff --git a/Code/Source/svFSILS/ns_solver.cpp b/Code/Source/svFSILS/ns_solver.cpp index 398a933..b05f3ee 100644 --- a/Code/Source/svFSILS/ns_solver.cpp +++ b/Code/Source/svFSILS/ns_solver.cpp @@ -270,19 +270,19 @@ void ns_solver(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_l // P = D*U // - auto P_col = P.col(i); - spar_mul::fsils_spar_mul_vs(lhs, lhs.rowPtr, lhs.colPtr, nsd, mD, U.slice(i), P_col); - P.set_col(i, P_col); + auto P_col = P.rcol(i); + spar_mul::fsils_spar_mul_vs(lhs, lhs.rowPtr, lhs.colPtr, nsd, mD, U.rslice(i), P_col); + //P.set_col(i, P_col); // P = Rc - P // - P.set_col(i, Rc - P.col(i)); + P.set_col(i, Rc - P_col); // P = [L + G^t*G]^-1*P // - P_col = P.col(i); + P_col = P.rcol(i); cgrad::schur(lhs, ls.CG, nsd, Gt, mG, mL, P_col); - P.set_col(i, P_col); + //P.set_col(i, P_col); // MU1 = G*P // @@ -290,10 +290,10 @@ void ns_solver(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_l dmsg << "i: " << i+1; dmsg << "iB: " << iB+1; #endif - P_col = P.col(i); - auto MU_iB = MU.slice(iB); + P_col = P.rcol(i); + auto MU_iB = MU.rslice(iB); spar_mul::fsils_spar_mul_sv(lhs, lhs.rowPtr, lhs.colPtr, nsd, mG, P_col, MU_iB); - MU.set_slice(iB, MU_iB); + //MU.set_slice(iB, MU_iB); // MU2 = Rm - G*P // @@ -302,29 +302,29 @@ void ns_solver(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_l // U = inv(K) * [Rm - G*P] // lhs.debug_active = true; - auto U_i = U.slice(i); + auto U_i = U.rslice(i); gmres::gmres(lhs, ls.GM, nsd, mK, MU.slice(iBB), U_i); - U.set_slice(i, U_i); + //U.set_slice(i, U_i); // MU2 = K*U // - auto MU_iBB = MU.slice(iBB); - spar_mul::fsils_spar_mul_vv(lhs, lhs.rowPtr, lhs.colPtr, nsd, mK, U.slice(i), MU_iBB); - MU.set_slice(iBB, MU_iBB); + auto MU_iBB = MU.rslice(iBB); + spar_mul::fsils_spar_mul_vv(lhs, lhs.rowPtr, lhs.colPtr, nsd, mK, U.rslice(i), MU_iBB); + //MU.set_slice(iBB, MU_iBB); - add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_ADD, nsd, U.slice(i), MU_iBB); - MU.set_slice(iBB, MU_iBB); + add_bc_mul::add_bc_mul(lhs, BcopType::BCOP_TYPE_ADD, nsd, U.rslice(i), MU_iBB); + //MU.set_slice(iBB, MU_iBB); // MP1 = L*P // - auto MP_iB = MP.col(iB); - spar_mul::fsils_spar_mul_ss(lhs, lhs.rowPtr, lhs.colPtr, mL, P.col(i), MP_iB); - MP.set_col(iB, MP_iB); + auto MP_iB = MP.rcol(iB); + spar_mul::fsils_spar_mul_ss(lhs, lhs.rowPtr, lhs.colPtr, mL, P.rcol(i), MP_iB); + //MP.set_col(iB, MP_iB); // MP2 = D*U - auto MP_iBB = MP.col(iBB); - spar_mul::fsils_spar_mul_vs(lhs, lhs.rowPtr, lhs.colPtr, nsd, mD, U.slice(i), MP_iBB); - MP.set_col(iBB, MP_iBB); + auto MP_iBB = MP.rcol(iBB); + spar_mul::fsils_spar_mul_vs(lhs, lhs.rowPtr, lhs.colPtr, nsd, mD, U.rslice(i), MP_iBB); + //MP.set_col(iBB, MP_iBB); int c = 0;