Commit f9ea0c46 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Exploit symmetry of strain matrices

Measurements indicate speedup around 15% for a 3d brittle-fracture
simulation.
parent ed739722
Pipeline #8667 passed with stage
in 3 minutes and 52 seconds
......@@ -196,9 +196,9 @@ public:
// mechanical part of gradient psi,{eps} in physical space
Range gradient(0.0);
for (std::size_t i=0; i<dim; i++)
for (std::size_t j=0; j<dim; j++)
for (std::size_t j=0; j<=i; j++)
for (std::size_t a=0; a<dim; a++)
gradient[linearIndex(j,i)] += gradientSpec[a]*eigenV[a][i]*eigenV[a][j];
gradient[linearIndex(j,i)] += ((i==j) ? 1 : 2) * gradientSpec[a]*eigenV[a][i]*eigenV[a][j];
return gradient;
}
......@@ -259,9 +259,9 @@ public:
Range result(0.0);
for (std::size_t i=0;i < dim; ++i)
for (std::size_t j=0;j < dim; ++j)
for (std::size_t j=0;j <= i; ++j)
for (std::size_t k=0;k < dim; ++k)
for (std::size_t l=0;l < dim; ++l)
for (std::size_t l=0;l <= k; ++l)
{
// The result matrix is symmetric: Compute only the lower triangular part
if (linearIndex(i,j) < linearIndex(k,l))
......@@ -290,6 +290,12 @@ public:
}
}
// Off-diagonal entries count twice: They come in identical pairs,
// but we compute only one of them.
const auto ijFactor = (i==j) ? 1 : 2;
const auto klFactor = (k==l) ? 1 : 2;
value *= ijFactor * klFactor;
// Add new value
result[linearIndex(i,j)][linearIndex(k,l)] += value;
if (linearIndex(i,j) !=linearIndex(k,l)) // Off-diagonal entry?
......@@ -354,9 +360,9 @@ public:
// psi,{eps eps} in physical space
for (std::size_t i=0;i < dim; ++i)
for (std::size_t j=0;j < dim; ++j)
for (std::size_t j=0;j <= i; ++j)
for (std::size_t k=0;k < dim; ++k)
for (std::size_t l=0;l < dim; ++l)
for (std::size_t l=0;l <= k; ++l)
{
// The result matrix is symmetric: Compute only the lower triangular part
if (linearIndex(i,j) < linearIndex(k,l))
......@@ -385,6 +391,12 @@ public:
}
}
// Off-diagonal entries count twice: They come in identical pairs,
// but we compute only one of them.
auto ijFactor = (i==j) ? 1 : 2;
auto klFactor = (k==l) ? 1 : 2;
value *= ijFactor * klFactor;
// Add new value
result[linearIndex(i,j)][linearIndex(k,l)] += value;
if (linearIndex(i,j) !=linearIndex(k,l)) // Off-diagonal entry?
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment