diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp index 2c73633..b6e5217 100644 --- a/mgrid_2/getmatrix.cpp +++ b/mgrid_2/getmatrix.cpp @@ -967,15 +967,15 @@ void CalcElem_MasseSpecific(int const ial[3], double const xc[], double const c, //x32 = xc[i2 + 0] - xc[i3 + 0], y32 = xc[i2 + 1] - xc[i3 + 1]; const double jac = fabs(x21 * y13 - x13 * y21); - ske[0][0] += c * jac / 12.0; - ske[0][1] += c * jac / 24.0; - ske[0][2] += c * jac / 24.0; - ske[1][0] += c * jac / 24.0; - ske[1][1] += c * jac / 12.0; - ske[1][2] += c * jac / 24.0; - ske[2][0] += c * jac / 24.0; - ske[2][1] += c * jac / 24.0; - ske[2][2] += c * jac / 12.0; + ske[0][0] = c * jac / 12.0; + ske[0][1] = c * jac / 24.0; + ske[0][2] = c * jac / 24.0; + ske[1][0] = c * jac / 24.0; + ske[1][1] = c * jac / 12.0; + ske[1][2] = c * jac / 24.0; + ske[2][0] = c * jac / 24.0; + ske[2][1] = c * jac / 24.0; + ske[2][2] = c * jac / 12.0; return; } diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp index 91eb1b3..28e27cd 100644 --- a/mgrid_2/main.cpp +++ b/mgrid_2/main.cpp @@ -31,7 +31,7 @@ int main(int argc, char **argv ) // ########################################## double dt = 1.0; // time step - int steps = 1; // number of time iterations + int steps = 20; // number of time iterations double u0_mug = 18.0; double u0_fluid = 80.0; @@ -52,7 +52,7 @@ int main(int argc, char **argv ) SK.CalculateLaplace_mult(fv); // stiffness matrix (+K) SK.AddMass_mult(fv, 1.0/dt); // add mass matrix (+M/dt) - SK.ApplyRobinBC_mult(fv, u_out); // apply Robin-bnd (+C = +F) + // SK.ApplyRobinBC_mult(fv, u_out); // apply Robin-bnd (+C = +F) // SK = M/dt + K + C = F // SK.Debug(); @@ -75,26 +75,29 @@ int main(int argc, char **argv ) mesh_c.Init_Solution_mult(uv, 0, [u0_mug](double x, double y) -> double { return u0_mug; }); // mug mesh_c.Init_Solution_mult(uv, 1, [u0_fluid](double x, double y) -> double { return u0_fluid; }); // fluid mesh_c.Init_Solution_mult(uv, 2, [u0_air](double x, double y) -> double { return u0_air; }); // air + - // TODO DINO auto t3 = system_clock::now(); // start timer for (int step = 0; step < steps; ++step) { vector G(Mdt.Nrows(), 0.0); Mdt.Mult(G, uv); // G = M/dt * u_{n} + + vector H = fv; for (size_t i = 0; i < Mdt.Nrows(); ++i) { - fv[i] += G[i]; // F + G + H[i] += G[i]; // H = F + G } - JacobiSolve(SK, fv, uv); // solve: (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n} + JacobiSolve(SK, H, uv); // solve: (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n} + // ----- SK ----- ------ H ------- } auto t4 = system_clock::now(); // stop timer auto duration = duration_cast(t4 - t3); // duration in microseconds double t_diff = static_cast(duration.count()) / 1e6; // overall duration in seconds - cout << "JacobiSolve: timing in sec. : " << t_diff << endl; + cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl; mesh_c.Visualize(uv); diff --git a/mgrid_2/uv.txt b/mgrid_2/uv.txt index 85fe8fd..543db25 100644 --- a/mgrid_2/uv.txt +++ b/mgrid_2/uv.txt @@ -710,263 +710,263 @@ 51 258 259 198 249 260 224 253 260 -18 -18 -18 -18 -18 -80 -80 -18 -18 -18 -18 -18 -80 -80 -80 -80 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -18 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -18 -18 -18 -80 -80 -80 -18 -80 -18 -18 -18 -18 -18 -18 -18 -80 -80 -80 -80 -80 -80 -80 -18 -18 -18 -18 -18 -80 -18 -18 -80 -80 -18 -80 -80 -18 -18 -18 -80 -80 -18 -18 -18 -18 -18 -18 -80 -18 -18 -18 -80 -80 -80 -80 -18 -18 -80 -80 -80 -80 -18 -18 -80 -18 -18 -80 -80 -18 -80 -80 -80 -18 -18 -18 -18 -18 -80 -18 -80 -80 -80 -80 -80 -80 -18 -80 -18 -80 -80 -18 -18 -80 -80 -18 -80 -18 -18 -18 -18 -80 -18 -80 -18 -18 -18 -80 -80 -80 -80 -80 -18 -80 -18 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 -80 +56.5253 +56.3976 +18.0028 +18.0036 +41.3884 +59.2188 +59.481 +41.4993 +18.0043 +18.0033 +64.8405 +64.8559 +65.6073 +65.4396 +66.6363 +66.5438 +58.7974 +62.486 +62.7851 +61.9604 +58.4362 +58.9828 +65.3201 +66.1555 +64.6712 +64.6552 +65.1662 +65.5044 +65.4384 +18.0044 +18.0061 +18.0064 +17.9752 +18.0019 +19.6357 +25.1439 +32.4962 +48.0285 +56.1788 +63.1956 +65.9046 +62.0992 +65.2948 +64.9761 +65.1454 +62.8772 +65.8154 +67.781 +66.1349 +64.9959 +65.1175 +32.5673 +25.1803 +19.6487 +18.0042 +17.9768 +18.0076 +18.007 +18.0051 +18.0036 +18.0048 +18.006 +17.9871 +17.9179 +19.1685 +29.0629 +45.9221 +58.1109 +63.1187 +48.7956 +48.6474 +42.9161 +42.8779 +37.6987 +27.3217 +26.5918 +33.9646 +42.2497 +48.6139 +48.5933 +18.2339 +18.5746 +19.0268 +19.5599 +20.0597 +20.3387 +20.3654 +20.1263 +19.7363 +19.3129 +18.8616 +18.4843 +18.1976 +65.4946 +65.6116 +65.1553 +64.4327 +63.993 +66.4169 +66.2598 +59.2317 +63.1182 +58.0294 +45.8038 +29.0096 +19.1607 +17.9168 +17.986 +18.0052 +18.0042 +18.0031 +65.4314 +65.6788 +66.8552 +67.0153 +64.3896 +65.5502 +65.5754 +66.1087 +66.7257 +65.8483 +63.2582 +56.2835 +48.1449 +66.5784 +66.1519 +65.7913 +65.7752 +79.988 +22.6127 +24.0142 +19.7519 +81.4422 +80.989 +84.5671 +29.8613 +84.1789 +18.8819 +64.8316 +18.2947 +64.8163 +64.4568 +36.8334 +36.7422 +61.6727 +61.6764 +79.6179 +81.3844 +81.3569 +79.8651 +79.9072 +21.3818 +26.231 +22.3415 +28.2768 +21.1454 +80.799 +19.7983 +32.5731 +80.107 +80.5701 +25.5144 +66.7723 +68.441 +34.6158 +19.0523 +18.7985 +77.8373 +78.3352 +19.3092 +18.408 +64.696 +36.9817 +64.9483 +36.8905 +80.3483 +20.7322 +25.4617 +20.8238 +79.7957 +82.348 +82.3798 +79.2733 +31.7661 +21.6786 +79.7477 +75.9929 +81.9226 +73.538 +26.4093 +22.2508 +80.6383 +35.9364 +19.7953 +78.7616 +79.9034 +37.3081 +77.7677 +67.3328 +66.7926 +18.4607 +18.6238 +18.401 +18.8327 +64.3039 +79.8238 +20.7796 +80.0313 +79.8331 +80.2912 +80.1094 +81.4642 +80.453 +27.9032 +79.6193 +23.8404 +72.2419 +80.7603 +20.0339 +20.6096 +79.5425 +80.5261 +23.2727 +80.5126 +19.6855 +27.5709 +19.5357 +36.0337 +80.4916 +24.3447 +65.1814 +21.3227 +35.8812 +18.7452 +76.7136 +74.4788 +74.3452 +73.2167 +81.3952 +21.3436 +80.0416 +20.7984 +79.9301 +78.9373 +79.8146 +79.7826 +80.244 +79.8804 +80.3524 +80.3655 +78.7719 +80.304 +77.8566 +75.9527 +80.0566 +80.3426 +80.1979 +76.6955 +75.2345 +78.7903