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

mass matrix

Browse source
This commit is contained in:
dino.celebic 2026-01-24 15:01:54 +01:00
commit 02072fde90
3 changed files with 88 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

69
mgrid_2/getmatrix.cpp
View file

@ -456,30 +456,66 @@ double FEM_Matrix::ThermalConductivity(const int subdomain)
double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
{
double lambda = 0.0;
double c = 0.0;
switch (subdomain)
{
// ceramic mug
case 1:
lambda = 2.0 * 1e6;
c = 2.0 * 1e6;
break;
// water
case 2:
lambda = 4.184 * 1e6;
c = 4.184 * 1e6;
break;
// air
case 3:
lambda = 1.2 * 1e3;
c = 1.2 * 1e3;
break;
default:
lambda = 1.0;
c = 1.0;
break;
}
return lambda;
return c;
}
void FEM_Matrix::AddMass_mult(vector<double> &f)
{
cout << "\n############ FEM_Matrix::AddMass_mult ";
double tstart = omp_get_wtime(); // OpenMP
assert(_mesh.NdofsElement() == 3); // only for triangular, linear elements
//cout << _nnz << " vs. " << _id[_nrows] << " " << _nrows<< endl;
assert(_nnz == _id[_nrows]);
double ske[3][3], fe[3];
// Loop over all elements
auto const nelem = _mesh.Nelems();
auto const &ia = _mesh.GetConnectivity();
auto const &xc = _mesh.GetCoords();
const vector<int> sd_vec = _mesh.ElementSubdomains;
#pragma omp parallel for private(ske,fe)
for (int i = 0; i < nelem; ++i) {
auto subdomain = sd_vec[i];
double c = VolumetricHeatCapacity(subdomain);
//cout << subdomain << endl;
CalcElem_MasseSpecific(ia.data() + 3 * i, xc.data(), c, ske);
//AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
AddElem_3(ia.data() + 3 * i, ske, fe, f);
}
double duration = omp_get_wtime() - tstart; // OpenMP
cout << "finished in " << duration << " sec. ########\n"; // ToDo: change to systemclock
//Debug();
return;
}
void FEM_Matrix::CalculateLaplace(vector<double> &f)
@ -932,6 +968,27 @@ void CalcElem_Masse(int const ial[3], double const xc[], double ske[3][3])
return;
}
void CalcElem_MasseSpecific(int const ial[3], double const xc[], double const c, double ske[3][3])
{
const int i1 = 2 * ial[0], i2 = 2 * ial[1], i3 = 2 * ial[2];
const double x13 = xc[i3 + 0] - xc[i1 + 0], y13 = xc[i3 + 1] - xc[i1 + 1],
x21 = xc[i1 + 0] - xc[i2 + 0], y21 = xc[i1 + 1] - xc[i2 + 1];
//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;
return;
}
// #####################################################################
BisectInterpolation::BisectInterpolation()