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

changed: simplify some CompatibleOperators

Browse Source 
use add and outer_product instead of hand-rolled loops
This commit is contained in:
Arne Morten Kvarving
2018-11-30 13:20:07 +01:00
parent dc9b010d9f
commit 417525cd3b
1 changed files with 12 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
Apps/Common/CompatibleOperators.C
View File

@@ -22,10 +22,9 @@ 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 i = 1; i <= fe.basis(n).size(); ++i)
for (size_t j = 1; j <= fe.basis(n).size(); ++j)
for (size_t k = 1; k <= nsd; ++k)
EM[n](i,j) += scale*AC[k]*fe.grad(n)(j,k)*fe.basis(n)(i)*fe.detJxW;
for (size_t k = 1; k <= nsd; ++k)
EM[n].outer_product(fe.basis(n), fe.grad(n).getColumn(k),
true,scale*AC[k]*fe.detJxW);
}
@@ -59,9 +58,8 @@ void CompatibleOperators::Weak::Gradient(std::vector<Matrix>& EM,
{
size_t nsd = fe.grad(1).cols();
for (size_t n = 1; n <= nsd; ++n)
for (size_t i=1; i <= fe.basis(n).size(); ++i)
for (size_t j=1; j <= fe.basis(nsd+1).size(); ++j)
EM[10+5*(n-1)](i,j) += -scale*fe.grad(n)(i,n)*fe.basis(nsd+1)(j)*fe.detJxW;
EM[10+5*(n-1)].outer_product(fe.grad(n).getColumn(n), fe.basis(nsd+1),
true, -scale*fe.detJxW);
}
@@ -69,8 +67,7 @@ void CompatibleOperators::Weak::Gradient(Vectors& EV, const FiniteElement& fe,
double scale)
{
for (size_t n = 1; n <= fe.grad(n).cols(); ++n)
for (size_t i = 1; i <= fe.basis(n).size(); ++i)
EV[n](i) += scale*fe.grad(n)(i,n)*fe.detJxW;
EV[n].add(fe.grad(n).getColumn(n), scale*fe.detJxW);
}
@@ -83,12 +80,11 @@ void CompatibleOperators::Weak::Laplacian(std::vector<Matrix>& EM,
EqualOrderOperators::Weak::Laplacian(EM[n], fe, scale, false, n);
for (size_t m = 1; m <= nsd && stress; m++)
for (size_t i = 1; i <= fe.basis(m).size(); ++i)
for (size_t n = m; n <= nsd; n++)
for (size_t j = 1; j <= fe.basis(n).size(); ++j) {
int idx = m == n ? m : (m == 1 ? 6+n-m : 13);
EM[idx](i,j) += scale* fe.grad(n)(j,m) * fe.grad(m)(i,n) * fe.detJxW;
}
for (size_t n = m; n <= nsd; n++) {
int idx = m == n ? m : (m == 1 ? 6+n-m : 13);
EM[idx].outer_product(fe.grad(m).getColumn(n),
fe.grad(n).getColumn(m), true, scale*fe.detJxW);
}
}
@@ -127,8 +123,7 @@ void CompatibleOperators::Residual::Convection(Vectors& EV, const FiniteElement&
double conv = 0.0;
for (size_t k = 1; k<= nsd; ++k)
conv += U[k-1]*dUdX(n, k);
for (size_t i=1; i <= fe.basis(n).size(); ++i)
EV[n](i) -= scale * conv * fe.basis(n)(i) * fe.detJxW;
EV[n].add(fe.basis(n), -scale*conv*fe.detJxW);
}
}