OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add convection/advection variants for div-compatible

Browse Source
This commit is contained in:
Arne Morten Kvarving 2022-02-22 12:28:15 +01:00
parent b467fa1ca2
commit bd2b2b88f6
1 changed files with 62 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
Apps/Common/CompatibleOperators.C
View File

@ -24,8 +24,18 @@ void CompatibleOperators::Weak::Advection(std::vector<Matrix>& EM,
size_t nsd = fe.grad(1).cols();
for (size_t n = 1; n <= nsd; ++n)
for (size_t k = 1; k <= nsd; ++k)
if (cnvForm == WeakOperators::CONVECTIVE)
EM[n].outer_product(fe.basis(n), fe.grad(n).getColumn(k),
true,scale*AC[k-1]*fe.detJxW);
else if (cnvForm == WeakOperators::CONSERVATIVE)
EM[n].outer_product(fe.grad(n).getColumn(k),
fe.basis(n), true, -scale*AC[k-1]*fe.detJxW);
else if (cnvForm == WeakOperators::SKEWSYMMETRIC) {
EM[n].outer_product(fe.basis(n), fe.grad(n).getColumn(k),
true,0.5*scale*AC[k-1]*fe.detJxW);
EM[n].outer_product(fe.grad(n).getColumn(k),
fe.basis(n), true, -0.5*scale*AC[k-1]*fe.detJxW);
}
}
@ -40,16 +50,31 @@ void CompatibleOperators::Weak::Convection(std::vector<Matrix>& EM,
{12, 2, 13},
{17, 18, 3}};
size_t nsd = fe.grad(1).cols();
for (size_t m=1; m <= nsd; ++m)
for (size_t i=1; i <= fe.basis(m).size(); ++i)
for (size_t n = 1; n <= nsd; ++n)
for (size_t j=1; j <= fe.basis(n).size(); ++j) {
double conv = fe.basis(n)(j) * dUdX(m,n);
if (m == n)
for (size_t k = 1; k <= nsd; ++k)
conv += U[k-1] * fe.grad(n)(j,k);
EM[vidx[m-1][n-1]](i,j) += scale * conv * fe.basis(m)(i) * fe.detJxW;
}
for (size_t k = 1; k <= nsd; ++k)
for (size_t l = 1; l <= nsd; ++l) {
switch (form) {
case WeakOperators::CONVECTIVE:
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.basis(l), true, scale*fe.detJxW*dUdX(k,l));
if(k == l)
for (size_t m = 1; m <= nsd; ++m)
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.grad(l).getColumn(m), true, scale*fe.detJxW*U[m-1]);
break;
case WeakOperators::CONSERVATIVE:
EM[vidx[k-1][l-1]].outer_product(fe.grad(k).getColumn(l), fe.basis(l), true, -U[k-1]*scale*fe.detJxW);
if (k == l)
for (size_t m = 1; m <= nsd; ++m)
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.grad(l).getColumn(m), true, -U[m-1]*scale*fe.detJxW);
break;
case WeakOperators::SKEWSYMMETRIC:
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.basis(l), true, 0.5*scale*fe.detJxW*dUdX(k,l));
EM[vidx[k-1][l-1]].outer_product(fe.grad(k).getColumn(l), fe.basis(l), true, -0.5*U[k-1]*scale*fe.detJxW);
if (k == l)
for (size_t m = 1; m <= nsd; ++m) {
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.grad(l).getColumn(m), true, 0.5*scale*fe.detJxW*U[m-1]);
EM[vidx[k-1][l-1]].outer_product(fe.basis(k), fe.grad(l).getColumn(m), true, -0.5*U[m-1]*scale*fe.detJxW);
}
}
}
}
@ -120,11 +145,34 @@ void CompatibleOperators::Residual::Convection(Vectors& EV, const FiniteElement&
WeakOperators::ConvectionForm form)
{
size_t nsd = fe.grad(1).cols();
for (size_t n=1; n <= nsd; ++n) {
for (size_t k = 1; k <= nsd; ++k) {
double conv = 0.0;
for (size_t k = 1; k<= nsd; ++k)
conv += U[k-1]*dUdX(n, k);
EV[n].add(fe.basis(n), -scale*conv*fe.detJxW);
switch (form) {
case WeakOperators::CONVECTIVE:
for (size_t l = 1; l <= nsd; ++l)
conv += -UC[l-1]*dUdX(k,l);
EV[k].add(fe.basis(k), scale*conv*fe.detJxW);
break;
case WeakOperators::CONSERVATIVE:
for (size_t i = 1; i <= fe.basis(k).size(); ++i) {
conv = 0.0;
for (size_t l = 1; l <= nsd; ++l)
conv += U[k-1]*UC[l-1]*fe.grad(k)(i,l);
EV[k](i) += conv*scale*fe.detJxW;
}
break;
case WeakOperators::SKEWSYMMETRIC:
for (size_t l = 1; l <= nsd; ++l)
conv += -UC[l-1]*dUdX(k,l);
EV[k].add(fe.basis(k), 0.5*scale*conv*fe.detJxW);
for (size_t i = 1; i <= fe.basis(k).size(); ++i) {
conv = 0.0;
for (size_t l = 1; l <= nsd; ++l)
conv += U[k-1]*UC[l-1]*fe.grad(k)(i,l);
EV[k](i) += 0.5*scale*conv*fe.detJxW;
}
break;
}
}
}