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

changed: optimize ResidualOperators::Laplacian

Browse Source 
and remove templating for the EqualOrder one
This commit is contained in:
Arne Morten Kvarving 2023-10-25 11:30:08 +02:00
parent b02de8d87b
commit c4ea583708
3 changed files with 31 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Apps/Common/CompatibleOperators.C
View File

@ -183,11 +183,13 @@ void CompatibleOperators::Residual::Laplacian(Vectors& EV,
double scale, bool stress)
{
size_t nsd = fe.grad(1).cols();
for (size_t k = 1; k <= nsd; ++k)
for (size_t i = 1; i <= fe.basis(k).size(); ++i) {
double diff = 0.0;
for (size_t m = 1; m <= nsd; ++m)
diff += fe.grad(k)(i,m)*dUdX(k,m);
EV[k](i) += scale*diff*fe.detJxW;
}
auto dUdXT = dUdX;
dUdXT.transpose();
for (size_t k = 1; k <= nsd; ++k) {
Vector diff;
fe.grad(k).multiply(Vector(dUdXT[k-1].ptr(), nsd), diff);
if (stress)
fe.grad(k).multiply(Vector(dUdX[k-1].ptr(), nsd), diff, false, 1);
EV[k].add(diff, scale*fe.detJxW);
}
}

18
Apps/Common/EqualOrderOperators.C
View File

@ -302,3 +302,21 @@ void EqualOrderOperators::Residual::Laplacian (Vector& EV,
size_t nsd = fe.grad(basis).cols();
fe.grad(basis).multiply(Vector(dUdX.ptr(),nsd),EV,scale*fe.detJxW,1.0);
}
void EqualOrderOperators::Residual::Laplacian (Vector& EV,
const FiniteElement& fe,
const Tensor& dUdX,
double scale,
bool stress, int basis)
{
size_t nsd = fe.grad(1).cols();
auto dUdXT = dUdX;
dUdXT.transpose();
for (size_t k = 1; k <= nsd; ++k) {
Vector diff;
fe.grad(basis).multiply(Vector(dUdXT[k-1].ptr(), nsd), diff);
if (stress)
fe.grad(basis).multiply(Vector(dUdX[k-1].ptr(), nsd), diff, false, 1);
EV.add(diff, scale*fe.detJxW, 0, 1, k-1, nsd);
}
}

32
Apps/Common/EqualOrderOperators.h
View File

@ -191,8 +191,8 @@ public:
//! \param[in] scale Scaling factor for contribution
//! \param[in] basis Basis to use
static void Laplacian(Vector& EV, const FiniteElement& fe,
const Vec3& dUdX, double scale=1.0,
int basis=1);
const Vec3& dUdX, double scale = 1.0,
int basis = 1);
//! \brief Compute a laplacian term in a residual vector.
//! \param EV The element vector to add contribution to
@ -201,33 +201,9 @@ public:
//! \param[in] scale Scaling factor for contribution
//! \param[in] stress Whether to add extra stress formulation terms
//! \param[in] basis Basis to use
template<class T>
static void Laplacian(Vector& EV, const FiniteElement& fe,
const T& dUdX, double scale=1.0,
bool stress=false, int basis=1)
{
size_t cmp = EV.size() / fe.basis(basis).size();
for (size_t i = 1;i <= fe.basis(basis).size();i++) {
for (size_t k = 1;k <= cmp;k++) {
double diff = 0.0;
for (size_t l = 1;l <= fe.grad(basis).cols();l++)
diff += dUdX(k,l)*fe.grad(basis)(i,l);
diff *= scale*fe.detJxW;
// Add residual to rhs of momentum equation
EV((i-1)*cmp + k) += diff;
}
}
// Use stress formulation
if (stress) {
for (size_t i = 1;i <= fe.basis(basis).size();i++)
for (size_t k = 1;k <= cmp;k++)
for (size_t l = 1;l <= cmp;l++)
// Diffusion
EV((i-1)*cmp + k) += scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW;
}
}
const Tensor& dUdX, double scale = 1.0,
bool stress = false, int basis = 1);
};
};