From d8ecb0bc06b86a3d5ddbabd7e7fdcd05b7edd355 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 5 Apr 2024 16:18:24 +0200
Subject: [PATCH] Use simpler matrix types for element Hesse matrices
Non-composite settings use Dune::Matrix<double>,
composite settings use a FieldMatrix where each entry
is Dune::Matrix<double>.
This is not quite as simple as the dune-functions convention,
which is to use Dune::Matrix in all cases. But it is still
simpler than before.
---
dune/gfe/assemblers/geodesicfeassembler.hh | 7 +++-
.../localgeodesicfeadolcstiffness.hh | 29 ++++++++-------
.../assemblers/localgeodesicfefdstiffness.hh | 6 +--
.../assemblers/localgeodesicfestiffness.hh | 14 ++-----
dune/gfe/assemblers/mixedgfeassembler.hh | 24 +++++++-----
test/localadolcstiffnesstest.cc | 37 ++++++++-----------
6 files changed, 58 insertions(+), 59 deletions(-)
diff --git a/dune/gfe/assemblers/geodesicfeassembler.hh b/dune/gfe/assemblers/geodesicfeassembler.hh
index 070c4164..b24be867 100644
--- a/dune/gfe/assemblers/geodesicfeassembler.hh
+++ b/dune/gfe/assemblers/geodesicfeassembler.hh
@@ -168,7 +168,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
localSolution[i] = sol[localView.index(i)];
std::vector<double> localGradient(numOfBaseFct*blocksize);
- Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > localHessian(numOfBaseFct,numOfBaseFct);
+ Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
@@ -181,7 +181,10 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
for (int j=0; j<numOfBaseFct; j++ ) {
auto col = localView.index(j);
- hessian[row][col] += localHessian[i][j];
+
+ for (std::size_t ii=0; ii<blocksize; ii++)
+ for (std::size_t jj=0; jj<blocksize; jj++)
+ hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 9220b8d1..594c0c50 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -390,7 +390,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
typedef typename TargetSpace::TangentVector TangentVector;
- localHessian.setSize(nDofs,nDofs);
+ localHessian.setSize(nDofs*blocksize,nDofs*blocksize);
for (size_t col=0; col<nDofs; col++) {
@@ -404,7 +404,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize; subRow++)
- localHessian[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[row*blocksize+subRow][col*blocksize+subCol] = semiEmbeddedProduct[subRow];
}
}
@@ -431,7 +431,8 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
TangentVector tmp2;
orthonormalFrame[row].mv(tmp1,tmp2);
- localHessian[row][row][subRow] += tmp2;
+ for (size_t subCol=0; subCol<blocksize; subCol++)
+ localHessian[row*blocksize+subRow][row*blocksize+subCol] += tmp2[subCol];
}
}
@@ -698,7 +699,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
using namespace Dune::Indices;
- localHessian[_0][_0].setSize(nDofs0,nDofs0);
+ localHessian[_0][_0].setSize(nDofs0*blocksize0, nDofs0*blocksize0);
for (size_t col=0; col<nDofs0; col++)
{
@@ -713,12 +714,12 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
- localHessian[_0][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_0][_0][row*blocksize0+subRow][col*blocksize0 + subCol] = semiEmbeddedProduct[subRow];
}
}
}
- localHessian[_0][_1].setSize(nDofs0,nDofs1);
+ localHessian[_0][_1].setSize(nDofs0*blocksize0, nDofs1*blocksize1);
for (size_t col=0; col<nDofs1; col++)
{
@@ -733,12 +734,12 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
- localHessian[_0][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_0][_1][row*blocksize0+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
}
}
}
- localHessian[_1][_0].setSize(nDofs1,nDofs0);
+ localHessian[_1][_0].setSize(nDofs1*blocksize1, nDofs0*blocksize0);
for (size_t col=0; col<nDofs0; col++)
{
@@ -752,12 +753,12 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian10[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
- localHessian[_1][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_1][_0][row*blocksize1+subRow][col*blocksize0+subCol] = semiEmbeddedProduct[subRow];
}
}
}
- localHessian[_1][_1].setSize(nDofs1,nDofs1);
+ localHessian[_1][_1].setSize(nDofs1*blocksize1, nDofs1*blocksize1);
for (size_t col=0; col<nDofs1; col++)
{
@@ -772,7 +773,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
- localHessian[_1][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_1][_1][row*blocksize1+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
}
}
}
@@ -802,7 +803,8 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typename TargetSpace0::TangentVector tmp2;
orthonormalFrame0[row].mv(tmp1,tmp2);
- localHessian[_0][_0][row][row][subRow] += tmp2;
+ for (size_t subCol=0; subCol<blocksize0; subCol++)
+ localHessian[_0][_0][row*blocksize0+subRow][row*blocksize0+subCol] += tmp2[subCol];
}
}
@@ -816,7 +818,8 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typename TargetSpace1::TangentVector tmp2;
orthonormalFrame1[row].mv(tmp1,tmp2);
- localHessian[_1][_1][row][row][subRow] += tmp2;
+ for (size_t subCol=0; subCol<blocksize1; subCol++)
+ localHessian[_1][_1][row*blocksize1+subRow][row*blocksize1+subCol] += tmp2[subCol];
}
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index b8c50b4c..630e9a1e 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -156,7 +156,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
size_t nDofs = localSolution.size();
// Clear assemble data
- localHessian.setSize(nDofs, nDofs);
+ localHessian.setSize(nDofs*blocksize, nDofs*blocksize);
localHessian = 0;
@@ -264,9 +264,9 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
#ifdef MULTIPRECISION
- localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo.template convert_to<double>();
+ localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo.template convert_to<double>();
#else
- localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo;
+ localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo;
#endif
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index b7848210..e1f936fa 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -28,10 +28,9 @@ namespace Dune::GFE
public:
// Type of the local Hessian
- using Hessian = Matrix<FieldMatrix<RT, blocksize, blocksize> >;
+ using Hessian = Matrix<double>;
- using Row = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize, blocksize> > >;
- using CompositeHessian = MultiTypeBlockMatrix<Row>;
+ using CompositeHessian = FieldMatrix<Matrix<double>,1,1>;
};
/** \brief A class exporting container types for sets local Hesse matrices
@@ -60,15 +59,10 @@ namespace Dune::GFE
public:
// Type of the local Hessian
- using Hessian = Matrix<FieldMatrix<RT, blocksize, blocksize> >;
+ using Hessian = Matrix<double>;
// Type of the local Hessian
- using Row0 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize0, blocksize0> >,
- Matrix<FieldMatrix<RT, blocksize0, blocksize1> > >;
- using Row1 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize1, blocksize0> >,
- Matrix<FieldMatrix<RT, blocksize1, blocksize1> > >;
-
- using CompositeHessian = MultiTypeBlockMatrix<Row0, Row1>;
+ using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
};
}
}
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index 4df0b50a..a3cef726 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -194,12 +194,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
std::vector<double> localGradient(nDofs0*blocksize0 + nDofs1*blocksize1);
- using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize0> >,
- Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > >;
- using Row1 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize1, blocksize0> >,
- Dune::Matrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > >;
-
- using HessianType = Dune::MultiTypeBlockMatrix<Row0, Row1>;
+ using HessianType = Dune::FieldMatrix<Dune::Matrix<double>,2,2>;
HessianType localHessian;
@@ -237,16 +232,25 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
auto col = localView.index(localIndexCol);
if (row[0]==0 and col[0]==0)
- hessian[_0][_0][row[1]][col[1]] += localHessian[_0][_0][i][j];
+ for (std::size_t ii=0; ii<blocksize0; ii++)
+ for (std::size_t jj=0; jj<blocksize0; jj++)
+ hessian[_0][_0][row[1]][col[1]][ii][jj] += localHessian[_0][_0][i*blocksize0+ii][j*blocksize0+jj];
if (row[0]==0 and col[0]==1)
- hessian[_0][_1][row[1]][col[1]] += localHessian[_0][_1][i][j-nDofs0];
+ for (std::size_t ii=0; ii<blocksize0; ii++)
+ for (std::size_t jj=0; jj<blocksize1; jj++)
+ hessian[_0][_1][row[1]][col[1]][ii][jj] += localHessian[_0][_1][i*blocksize0+ii][(j-nDofs0)*blocksize1+jj];
if (row[0]==1 and col[0]==0)
- hessian[_1][_0][row[1]][col[1]] += localHessian[_1][_0][i-nDofs0][j];
+ for (std::size_t ii=0; ii<blocksize1; ii++)
+ for (std::size_t jj=0; jj<blocksize0; jj++)
+ hessian[_1][_0][row[1]][col[1]][ii][jj] += localHessian[_1][_0][(i-nDofs0)*blocksize1+ii][j*blocksize0+jj];
if (row[0]==1 and col[0]==1)
- hessian[_1][_1][row[1]][col[1]] += localHessian[_1][_1][i-nDofs0][j-nDofs0];
+ for (std::size_t ii=0; ii<blocksize1; ii++)
+ for (std::size_t jj=0; jj<blocksize1; jj++)
+ hessian[_1][_1][row[1]][col[1]][ii][jj] += localHessian[_1][_1][(i-nDofs0)*blocksize1+ii][(j-nDofs0)*
+ blocksize1+jj];
}
// Add local gradient to global gradient
diff --git a/test/localadolcstiffnesstest.cc b/test/localadolcstiffnesstest.cc
index 5c5faef4..10ca626b 100644
--- a/test/localadolcstiffnesstest.cc
+++ b/test/localadolcstiffnesstest.cc
@@ -55,9 +55,8 @@ const int dim = 2;
using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
// Compare two matrices
-template <int N>
-void compareMatrices(const Matrix<FieldMatrix<double,N,N> >& matrixA, std::string nameA,
- const Matrix<FieldMatrix<double,N,N> >& matrixB, std::string nameB)
+void compareMatrices(const Matrix<double>& matrixA, std::string nameA,
+ const Matrix<double>& matrixB, std::string nameB)
{
double maxAbsDifference = -1;
double maxRelDifference = -1;
@@ -66,22 +65,18 @@ void compareMatrices(const Matrix<FieldMatrix<double,N,N> >& matrixA, std::strin
for (size_t j=0; j<matrixA.M(); j++ ) {
- for (size_t ii=0; ii<matrixA[i][j].N(); ii++)
- for (size_t jj=0; jj<matrixA[i][j].M(); jj++)
- {
- double valueA = matrixA[i][j][ii][jj];
- double valueB = matrixB[i][j][ii][jj];
-
- double absDifference = valueA - valueB;
- double relDifference = std::abs(absDifference) / std::abs(valueA);
- maxAbsDifference = std::max(maxAbsDifference, std::abs(absDifference));
- if (not std::isinf(relDifference))
- maxRelDifference = std::max(maxRelDifference, relDifference);
-
- if (relDifference > 1)
- std::cout << i << ", " << j << " " << ii << ", " << jj
- << ", " << nameA << ": " << valueA << ", " << nameB << ": " << valueB << std::endl;
- }
+ const double valueA = matrixA[i][j];
+ const double valueB = matrixB[i][j];
+
+ double absDifference = valueA - valueB;
+ double relDifference = std::abs(absDifference) / std::abs(valueA);
+ maxAbsDifference = std::max(maxAbsDifference, std::abs(absDifference));
+ if (not std::isinf(relDifference))
+ maxRelDifference = std::max(maxRelDifference, relDifference);
+
+ if (relDifference > 1)
+ std::cout << i << ", " << j << " , "
+ << nameA << ": " << valueA << ", " << nameB << ": " << valueB << std::endl;
}
}
@@ -213,8 +208,8 @@ int main (int argc, char *argv[]) try
std::vector<double> localRiemannianADGradient(numOfBaseFct*blocksize);
std::vector<double> localRiemannianFDGradient(numOfBaseFct*blocksize);
- Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianADHessian;
- Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianFDHessian;
+ Matrix<double> localRiemannianADHessian;
+ Matrix<double> localRiemannianFDHessian;
localGFEADOLCStiffness.assembleGradientAndHessian(localView,
localSolution,
--
GitLab