jakob.schratter/SciFEM-Project_CoffeeMugSimulation_Celebic-Schratter
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

fix bugs

Browse source
This commit is contained in:
dino.celebic 2026-01-25 13:14:37 +01:00
commit decfe1d710
3 changed files with 278 additions and 275 deletions
Download patch file Download diff file Expand all files Collapse all files

18
mgrid_2/getmatrix.cpp
View file

@ -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]; //x32 = xc[i2 + 0] - xc[i3 + 0], y32 = xc[i2 + 1] - xc[i3 + 1];
const double jac = fabs(x21 * y13 - x13 * y21); const double jac = fabs(x21 * y13 - x13 * y21);
ske[0][0] += c * jac / 12.0; ske[0][0] = c * jac / 12.0;
ske[0][1] += c * jac / 24.0; ske[0][1] = c * jac / 24.0;
ske[0][2] += c * jac / 24.0; ske[0][2] = c * jac / 24.0;
ske[1][0] += c * jac / 24.0; ske[1][0] = c * jac / 24.0;
ske[1][1] += c * jac / 12.0; ske[1][1] = c * jac / 12.0;
ske[1][2] += c * jac / 24.0; ske[1][2] = c * jac / 24.0;
ske[2][0] += c * jac / 24.0; ske[2][0] = c * jac / 24.0;
ske[2][1] += c * jac / 24.0; ske[2][1] = c * jac / 24.0;
ske[2][2] += c * jac / 12.0; ske[2][2] = c * jac / 12.0;
return; return;
} }

15
mgrid_2/main.cpp
View file

@ -31,7 +31,7 @@ int main(int argc, char **argv )
// ########################################## // ##########################################
double dt = 1.0; // time step 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_mug = 18.0;
double u0_fluid = 80.0; double u0_fluid = 80.0;
@ -52,7 +52,7 @@ int main(int argc, char **argv )
SK.CalculateLaplace_mult(fv); // stiffness matrix (+K) SK.CalculateLaplace_mult(fv); // stiffness matrix (+K)
SK.AddMass_mult(fv, 1.0/dt); // add mass matrix (+M/dt) 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 = M/dt + K + C = F
// SK.Debug(); // 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, 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, 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 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 auto t3 = system_clock::now(); // start timer
for (int step = 0; step < steps; ++step) for (int step = 0; step < steps; ++step)
{ {
vector<double> G(Mdt.Nrows(), 0.0); vector<double> G(Mdt.Nrows(), 0.0);
Mdt.Mult(G, uv); // G = M/dt * u_{n} Mdt.Mult(G, uv); // G = M/dt * u_{n}
vector<double> H = fv;
for (size_t i = 0; i < Mdt.Nrows(); ++i) 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 t4 = system_clock::now(); // stop timer
auto duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds auto duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds
double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds double t_diff = static_cast<double>(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); mesh_c.Visualize(uv);

520
mgrid_2/uv.txt
View file

@ -710,263 +710,263 @@
51 258 259 51 258 259
198 249 260 198 249 260
224 253 260 224 253 260
18 56.5253
18 56.3976
18 18.0028
18 18.0036
18 41.3884
80 59.2188
80 59.481
18 41.4993
18 18.0043
18 18.0033
18 64.8405
18 64.8559
80 65.6073
80 65.4396
80 66.6363
80 66.5438
18 58.7974
18 62.486
18 62.7851
18 61.9604
18 58.4362
18 58.9828
18 65.3201
18 66.1555
18 64.6712
18 64.6552
18 65.1662
18 65.5044
18 65.4384
18 18.0044
18 18.0061
18 18.0064
18 17.9752
18 18.0019
18 19.6357
18 25.1439
18 32.4962
80 48.0285
80 56.1788
80 63.1956
80 65.9046
80 62.0992
80 65.2948
80 64.9761
80 65.1454
80 62.8772
80 65.8154
80 67.781
80 66.1349
80 64.9959
80 65.1175
18 32.5673
18 25.1803
18 19.6487
18 18.0042
18 17.9768
18 18.0076
18 18.007
18 18.0051
18 18.0036
18 18.0048
18 18.006
18 17.9871
18 17.9179
18 19.1685
18 29.0629
18 45.9221
18 58.1109
18 63.1187
18 48.7956
18 48.6474
18 42.9161
18 42.8779
18 37.6987
18 27.3217
18 26.5918
18 33.9646
18 42.2497
18 48.6139
18 48.5933
18 18.2339
18 18.5746
18 19.0268
18 19.5599
18 20.0597
18 20.3387
18 20.3654
18 20.1263
18 19.7363
18 19.3129
18 18.8616
18 18.4843
18 18.1976
18 65.4946
18 65.6116
18 65.1553
18 64.4327
18 63.993
18 66.4169
18 66.2598
18 59.2317
18 63.1182
18 58.0294
18 45.8038
18 29.0096
18 19.1607
18 17.9168
18 17.986
18 18.0052
18 18.0042
18 18.0031
80 65.4314
80 65.6788
80 66.8552
80 67.0153
80 64.3896
80 65.5502
80 65.5754
80 66.1087
80 66.7257
80 65.8483
80 63.2582
80 56.2835
80 48.1449
80 66.5784
80 66.1519
80 65.7913
80 65.7752
80 79.988
18 22.6127
18 24.0142
18 19.7519
80 81.4422
80 80.989
80 84.5671
18 29.8613
80 84.1789
18 18.8819
18 64.8316
18 18.2947
18 64.8163
18 64.4568
18 36.8334
18 36.7422
80 61.6727
80 61.6764
80 79.6179
80 81.3844
80 81.3569
80 79.8651
80 79.9072
18 21.3818
18 26.231
18 22.3415
18 28.2768
18 21.1454
80 80.799
18 19.7983
18 32.5731
80 80.107
80 80.5701
18 25.5144
80 66.7723
80 68.441
18 34.6158
18 19.0523
18 18.7985
80 77.8373
80 78.3352
18 19.3092
18 18.408
18 64.696
18 36.9817
18 64.9483
18 36.8905
80 80.3483
18 20.7322
18 25.4617
18 20.8238
80 79.7957
80 82.348
80 82.3798
80 79.2733
18 31.7661
18 21.6786
80 79.7477
80 75.9929
80 81.9226
80 73.538
18 26.4093
18 22.2508
80 80.6383
18 35.9364
18 19.7953
80 78.7616
80 79.9034
18 37.3081
80 77.7677
80 67.3328
80 66.7926
18 18.4607
18 18.6238
18 18.401
18 18.8327
18 64.3039
80 79.8238
18 20.7796
80 80.0313
80 79.8331
80 80.2912
80 80.1094
80 81.4642
80 80.453
18 27.9032
80 79.6193
18 23.8404
80 72.2419
80 80.7603
18 20.0339
18 20.6096
80 79.5425
80 80.5261
18 23.2727
80 80.5126
18 19.6855
18 27.5709
18 19.5357
18 36.0337
80 80.4916
18 24.3447
80 65.1814
18 21.3227
18 35.8812
18 18.7452
80 76.7136
80 74.4788
80 74.3452
80 73.2167
80 81.3952
18 21.3436
80 80.0416
18 20.7984
80 79.9301
80 78.9373
80 79.8146
80 79.7826
80 80.244
80 79.8804
80 80.3524
80 80.3655
80 78.7719
80 80.304
80 77.8566
80 75.9527
80 80.0566
80 80.3426
80 80.1979
80 76.6955
80 75.2345
80 78.7903