diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index 9bbf4c0cc5535ab7cdde8a8e1ae5f5adbb75ab36..4b5382c21440a58caa205d0cf16513c06400e0cc 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -15,8 +15,6 @@
#include <dune/gfe/localenergy.hh>
#include <dune/gfe/localgeodesicfestiffness.hh>
-#define ADOLC_VECTOR_MODE
-
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class Basis, class TargetSpace>
@@ -42,8 +40,9 @@ public:
//! Dimension of the embedding space
constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
- LocalGeodesicFEADOLCStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy)
- : localEnergy_(energy)
+ LocalGeodesicFEADOLCStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy, bool adolcScalarMode = false)
+ : localEnergy_(energy),
+ adolcScalarMode_(adolcScalarMode)
{}
/** \brief Compute the energy at the current configuration */
@@ -67,7 +66,7 @@ public:
std::vector<typename TargetSpace::TangentVector>& localGradient);
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
-
+ const bool adolcScalarMode_;
};
@@ -251,62 +250,61 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
-#ifndef ADOLC_VECTOR_MODE
- std::vector<double> v(nDoubles);
- std::vector<double> w(nDoubles);
-
- std::fill(v.begin(), v.end(), 0.0);
-
- for (size_t i=0; i<nDofs; i++)
- for (int ii=0; ii<blocksize; ii++)
- {
- // Evaluate Hessian in the direction of each vector of the orthonormal frame
- for (size_t k=0; k<embeddedBlocksize; k++)
- v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
-
- int rc= 3;
- MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
- if (rc < 0)
- DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
-
- for (size_t j=0; j<nDoubles; j++)
- embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
-
- // Make v the null vector again
- std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
- }
-#else
- int n = nDoubles;
- int nDirections = nDofs * blocksize;
- double* tangent[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- tangent[i] = (double*)malloc(nDirections*sizeof(double));
-
- double* rawHessian[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
-
- for (int j=0; j<nDirections; j++)
- {
- for (int i=0; i<n; i++)
- tangent[i][j] = 0.0;
-
- for (int i=0; i<embeddedBlocksize; i++)
- tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
- }
- hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+ if (adolcScalarMode_) {
+ std::vector<double> v(nDoubles);
+ std::vector<double> w(nDoubles);
+
+ std::fill(v.begin(), v.end(), 0.0);
+
+ for (size_t i=0; i<nDofs; i++)
+ for (int ii=0; ii<blocksize; ii++)
+ {
+ // Evaluate Hessian in the direction of each vector of the orthonormal frame
+ for (size_t k=0; k<embeddedBlocksize; k++)
+ v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
+
+ int rc= 3;
+ MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
+ if (rc < 0)
+ DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+ for (size_t j=0; j<nDoubles; j++)
+ embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
+
+ // Make v the null vector again
+ std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
+ }
+ } else { //ADOL-C vector mode
+ int n = nDoubles;
+ int nDirections = nDofs * blocksize;
+ double* tangent[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ tangent[i] = (double*)malloc(nDirections*sizeof(double));
+
+ double* rawHessian[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
+
+ for (int j=0; j<nDirections; j++)
+ {
+ for (int i=0; i<n; i++)
+ tangent[i][j] = 0.0;
+
+ for (int i=0; i<embeddedBlocksize; i++)
+ tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
+ }
+ hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
- // Copy Hessian into Dune data type
- for(size_t i=0; i<nDoubles; i++)
- for (int j=0; j<nDirections; j++)
- embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
+ // Copy Hessian into Dune data type
+ for(size_t i=0; i<nDoubles; i++)
+ for (int j=0; j<nDirections; j++)
+ embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
- for(size_t i=0; i<nDoubles; i++) {
- free(rawHessian[i]);
- free(tangent[i]);
+ for(size_t i=0; i<nDoubles; i++) {
+ free(rawHessian[i]);
+ free(tangent[i]);
+ }
}
-#endif
-
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
diff --git a/dune/gfe/mixedlocalgfeadolcstiffness.hh b/dune/gfe/mixedlocalgfeadolcstiffness.hh
index 40afdad805870d4b037cee14fee930f8b761d229..92607e323d8cad340cb5b1180bf07cbb440bf270 100644
--- a/dune/gfe/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/mixedlocalgfeadolcstiffness.hh
@@ -14,8 +14,6 @@
#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
-#define ADOLC_VECTOR_MODE
-
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class Basis, class TargetSpace0, class TargetSpace1>
@@ -46,8 +44,9 @@ public:
constexpr static int embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension;
MixedLocalGFEADOLCStiffness(const MixedLocalGeodesicFEStiffness<Basis, ATargetSpace0,
- ATargetSpace1>* energy)
- : localEnergy_(energy)
+ ATargetSpace1>* energy, bool adolcScalarMode = false)
+ : localEnergy_(energy),
+ adolcScalarMode_(adolcScalarMode)
{}
/** \brief Compute the energy at the current configuration */
@@ -75,7 +74,7 @@ public:
std::vector<typename TargetSpace1::TangentVector>& localGradient1);
const MixedLocalGeodesicFEStiffness<Basis, ATargetSpace0, ATargetSpace1>* localEnergy_;
-
+ const bool adolcScalarMode_;
};
@@ -293,88 +292,118 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize0> > embeddedHessian10(nDofs1,nDofs0);
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize1> > embeddedHessian11(nDofs1,nDofs1);
-#ifndef ADOLC_VECTOR_MODE
-#error ADOL-C scalar mode not implemented
-#if 0
- std::vector<double> v(nDoubles);
- std::vector<double> w(nDoubles);
-
- std::fill(v.begin(), v.end(), 0.0);
-
- for (int i=0; i<nDofs; i++)
- for (int ii=0; ii<blocksize; ii++)
- {
- // Evaluate Hessian in the direction of each vector of the orthonormal frame
- for (size_t k=0; k<embeddedBlocksize; k++)
- v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
-
- int rc= 3;
- MINDEC(rc, hess_vec(1, nDoubles, xp.data(), v.data(), w.data()));
- if (rc < 0)
- DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
-
- for (int j=0; j<nDoubles; j++)
- embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
+ if(adolcScalarMode_) {
+ std::vector<double> v(nDoubles);
+ std::vector<double> w(nDoubles);
+
+ std::fill(v.begin(), v.end(), 0.0);
+
+ size_t nDoubles0 = nDofs0*embeddedBlocksize0; // nDoubles = nDoubles0 + nDoubles1
+ size_t nDoubles1 = nDofs1*embeddedBlocksize1;
+
+ std::fill(v.begin(), v.end(), 0.0);
+
+ for (size_t i=0; i<nDofs0 + nDofs1; i++) {
+
+ // Evaluate Hessian in the direction of each vector of the orthonormal frame
+ if (i < nDofs0) { //Upper half
+ auto i0 = i;
+ for (int ii0=0; ii0<blocksize0; ii0++) {
+ for (size_t k0=0; k0<embeddedBlocksize0; k0++) {
+ v[i0*embeddedBlocksize0 + k0] = orthonormalFrame0[i0][ii0][k0];
+ }
+ int rc= 3;
+ MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
+ if (rc < 0)
+ DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+ for (size_t j0=0; j0<nDoubles0; j0++) //Upper left
+ embeddedHessian00[i0][j0/embeddedBlocksize0][ii0][j0%embeddedBlocksize0] = w[j0];
+
+ for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
+ embeddedHessian01[i0][j1/embeddedBlocksize1][ii0][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
+
+ // Make v the null vector again
+ std::fill(&v[i0*embeddedBlocksize0], &v[(i0+1)*embeddedBlocksize0], 0.0);
+ }
+ } else { // i = nDofs0 ... nDofs0 + nDofs1 //lower half
+ auto i1 = i - nDofs0;
+ for (int ii1=0; ii1<blocksize1; ii1++) {
+ for (size_t k1=0; k1<embeddedBlocksize1; k1++) {
+ v[nDoubles0 + i1*embeddedBlocksize1 + k1] = orthonormalFrame1[i1][ii1][k1];
+ }
+ int rc= 3;
+ MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
+ if (rc < 0)
+ DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+ for (size_t j0=0; j0<nDoubles0; j0++) //Uppper left
+ embeddedHessian10[i1][j0/embeddedBlocksize0][ii1][j0%embeddedBlocksize0] = w[j0];
+
+ for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
+ embeddedHessian11[i1][j1/embeddedBlocksize1][ii1][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
+
+ // Make v the null vector again
+ std::fill(&v[nDoubles0 + i1*embeddedBlocksize1], &v[nDoubles0 + (i1+1)*embeddedBlocksize1], 0.0);
+ }
+ }
+ }
- // Make v the null vector again
- std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
- }
-#endif
-#else
- int n = nDoubles;
- int nDirections = nDofs0 * blocksize0 + nDofs1 * blocksize1;
- double* tangent[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- tangent[i] = (double*)malloc(nDirections*sizeof(double));
-
- double* rawHessian[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
-
- // Initialize directions field with zeros
- for (int j=0; j<nDirections; j++)
- for (int i=0; i<n; i++)
- tangent[i][j] = 0.0;
-
- for (size_t j=0; j<nDofs0*blocksize0; j++)
- for (int i=0; i<embeddedBlocksize0; i++)
- tangent[(j/blocksize0)*embeddedBlocksize0+i][j] = orthonormalFrame0[j/blocksize0][j%blocksize0][i];
-
- for (size_t j=0; j<nDofs1*blocksize1; j++)
- for (int i=0; i<embeddedBlocksize1; i++)
- tangent[nDofs0*embeddedBlocksize0 + (j/blocksize1)*embeddedBlocksize1+i][nDofs0*blocksize0 + j] = orthonormalFrame1[j/blocksize1][j%blocksize1][i];
-
- hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
-
- // Copy Hessian into Dune data type
- size_t offset0 = nDofs0*embeddedBlocksize0;
- size_t offset1 = nDofs0*blocksize0;
-
- // upper left block
- for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
- for (size_t j=0; j<nDofs0*blocksize0; j++)
- embeddedHessian00[j/blocksize0][i/embeddedBlocksize0][j%blocksize0][i%embeddedBlocksize0] = rawHessian[i][j];
-
- // upper right block
- for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
- for (size_t j=0; j<nDofs0*blocksize0; j++)
- embeddedHessian01[j/blocksize0][i/embeddedBlocksize1][j%blocksize0][i%embeddedBlocksize1] = rawHessian[offset0+i][j];
-
- // lower left block
- for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
- for (size_t j=0; j<nDofs1*blocksize1; j++)
- embeddedHessian10[j/blocksize1][i/embeddedBlocksize0][j%blocksize1][i%embeddedBlocksize0] = rawHessian[i][offset1+j];
-
- // lower right block
- for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
- for (size_t j=0; j<nDofs1*blocksize1; j++)
- embeddedHessian11[j/blocksize1][i/embeddedBlocksize1][j%blocksize1][i%embeddedBlocksize1] = rawHessian[offset0+i][offset1+j];
-
- for(size_t i=0; i<nDoubles; i++) {
- free(rawHessian[i]);
- free(tangent[i]);
+ } else { //ADOL-C vector mode}
+ int n = nDoubles;
+ int nDirections = nDofs0 * blocksize0 + nDofs1 * blocksize1;
+ double* tangent[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ tangent[i] = (double*)malloc(nDirections*sizeof(double));
+
+ double* rawHessian[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
+
+ // Initialize directions field with zeros
+ for (int j=0; j<nDirections; j++)
+ for (int i=0; i<n; i++)
+ tangent[i][j] = 0.0;
+
+ for (size_t j=0; j<nDofs0*blocksize0; j++)
+ for (int i=0; i<embeddedBlocksize0; i++)
+ tangent[(j/blocksize0)*embeddedBlocksize0+i][j] = orthonormalFrame0[j/blocksize0][j%blocksize0][i];
+
+ for (size_t j=0; j<nDofs1*blocksize1; j++)
+ for (int i=0; i<embeddedBlocksize1; i++)
+ tangent[nDofs0*embeddedBlocksize0 + (j/blocksize1)*embeddedBlocksize1+i][nDofs0*blocksize0 + j] = orthonormalFrame1[j/blocksize1][j%blocksize1][i];
+
+ hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+
+ // Copy Hessian into Dune data type
+ size_t offset0 = nDofs0*embeddedBlocksize0;
+ size_t offset1 = nDofs0*blocksize0;
+
+ // upper left block
+ for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
+ for (size_t j=0; j<nDofs0*blocksize0; j++)
+ embeddedHessian00[j/blocksize0][i/embeddedBlocksize0][j%blocksize0][i%embeddedBlocksize0] = rawHessian[i][j];
+
+ // upper right block
+ for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
+ for (size_t j=0; j<nDofs0*blocksize0; j++)
+ embeddedHessian01[j/blocksize0][i/embeddedBlocksize1][j%blocksize0][i%embeddedBlocksize1] = rawHessian[offset0+i][j];
+
+ // lower left block
+ for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
+ for (size_t j=0; j<nDofs1*blocksize1; j++)
+ embeddedHessian10[j/blocksize1][i/embeddedBlocksize0][j%blocksize1][i%embeddedBlocksize0] = rawHessian[i][offset1+j];
+
+ // lower right block
+ for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
+ for (size_t j=0; j<nDofs1*blocksize1; j++)
+ embeddedHessian11[j/blocksize1][i/embeddedBlocksize1][j%blocksize1][i%embeddedBlocksize1] = rawHessian[offset0+i][offset1+j];
+
+ for(size_t i=0; i<nDoubles; i++) {
+ free(rawHessian[i]);
+ free(tangent[i]);
+ }
}
-#endif
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index b4aa70abe91c47bd214e6c4a4ce9ea00cb538d67..cc0ea99152d42c3817531764e1f633f279721d74 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -139,6 +139,7 @@ int main (int argc, char *argv[]) try
const double mgTolerance = parameterSet.get<double>("mgTolerance");
const double baseTolerance = parameterSet.get<double>("baseTolerance");
const bool instrumented = parameterSet.get<bool>("instrumented");
+ const bool adolcScalarMode = parameterSet.get<bool>("adolcScalarMode", false);
std::string resultPath = parameterSet.get("resultPath", "");
// ///////////////////////////////////////
@@ -453,7 +454,7 @@ int main (int argc, char *argv[]) try
volumeLoad);
MixedLocalGFEADOLCStiffness<CompositeBasis,
RealTuple<double,3>,
- Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy);
+ Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy, adolcScalarMode);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,3>,
Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
@@ -551,7 +552,7 @@ int main (int argc, char *argv[]) try
&neumannBoundary,
neumannFunction,
volumeLoad);
- localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>>(&localCosseratEnergyPlanar);
+ localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>>(&localCosseratEnergyPlanar, adolcScalarMode);
GeodesicFEAssembler<DeformationFEBasis,TargetSpace> assembler(gridView, *localGFEADOLCStiffness);
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b292a33af2b6da6b5cdb459a61db0020f7a67a3a..8f62323c8462bace4cc680e2ff7b4da4bccb440d 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,5 +1,6 @@
set(TESTS
adolctest
+ adolctest-scalar-and-vector-mode
averagedistanceassemblertest
cosseratenergytest
cosseratrodenergytest
diff --git a/test/adolctest-scalar-and-vector-mode.cc b/test/adolctest-scalar-and-vector-mode.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2f4fffce180161186a8c908ef764d0dd940b61a4
--- /dev/null
+++ b/test/adolctest-scalar-and-vector-mode.cc
@@ -0,0 +1,222 @@
+#include <config.h>
+
+#include <array>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/tuplevector.hh>
+
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+// grid dimension
+const int gridDim = 2;
+
+// target dimension
+const int dim = 3;
+
+using namespace Dune;
+using namespace Indices;
+
+//differentiation method: ADOL-C
+using ValueType = adouble;
+
+//Types for the mixed space
+using DisplacementVector = std::vector<RealTuple<double,dim>>;
+using RotationVector = std::vector<Rotation<double,dim>>;
+using Vector = TupleVector<DisplacementVector, RotationVector>;
+const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
+using CorrectionTypeMixed = MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>>;
+
+using MatrixRow0 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dim,dim>>, BCRSMatrix<FieldMatrix<double,dim,dimCR>>>;
+using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>>, BCRSMatrix<FieldMatrix<double,dimCR,dimCR>>>;
+using MatrixTypeMixed = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
+
+//Types for the Non-mixed space
+using RBM = RigidBodyMotion<double, dim>;
+using RBMVector = std::vector<RBM>;
+const static int blocksize = RBM::TangentVector::dimension;
+using CorrectionType = BlockVector<FieldVector<double, blocksize> >;
+using MatrixType = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
+
+int main (int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create the grid
+ /////////////////////////////////////////////////////////////////////////
+ using GridType = UGGrid<gridDim>;
+ auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1,1}, {2,2});
+ grid->globalRefine(2);
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create a composite basis
+ /////////////////////////////////////////////////////////////////////////
+
+ using namespace Functions::BasisFactory;
+
+ auto compositeBasisMixed = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ lagrange<2>()
+ ),
+ power<dim>(
+ lagrange<1>()
+ )
+ ));
+
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ lagrange<2>()
+ ),
+ power<dim>(
+ lagrange<2>()
+ )
+ ));
+
+ using CompositeBasis = decltype(compositeBasis);
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create the energy functions with their parameters
+ /////////////////////////////////////////////////////////////////////////
+
+ //Surface-Cosserat-Energy-Parameters
+ ParameterTree parameters;
+ parameters["thickness"] = "1";
+ parameters["mu"] = "2.7191e+4";
+ parameters["lambda"] = "4.4364e+4";
+ parameters["mu_c"] = "0";
+ parameters["L_c"] = "0.01";
+ parameters["q"] = "2";
+ parameters["kappa"] = "1";
+
+ //Mixed space
+ CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergyMixed(parameters,
+ nullptr,
+ nullptr,
+ nullptr);
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedLocalGFEADOLCStiffnessVector(&cosseratEnergyMixed, false);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssemblerVector(compositeBasis, &mixedLocalGFEADOLCStiffnessVector);
+
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedLocalGFEADOLCStiffnessScalar(&cosseratEnergyMixed, true);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssemblerScalar(compositeBasis, &mixedLocalGFEADOLCStiffnessScalar);
+
+ //Non-mixed space
+ using DeformationFEBasis = Dune::Functions::LagrangeBasis<GridView,2>;
+
+ CosseratEnergyLocalStiffness<DeformationFEBasis, dim,adouble> cosseratEnergy(parameters,
+ nullptr,
+ nullptr,
+ nullptr);
+
+ LocalGeodesicFEADOLCStiffness<DeformationFEBasis,RBM> localGFEADOLCStiffnessVector(&cosseratEnergy, false);
+ GeodesicFEAssembler<DeformationFEBasis,RBM> assemblerVector(gridView, localGFEADOLCStiffnessVector);
+
+ LocalGeodesicFEADOLCStiffness<DeformationFEBasis,RBM> localGFEADOLCStiffnessScalar(&cosseratEnergy, true);
+ GeodesicFEAssembler<DeformationFEBasis,RBM> assemblerScalar(gridView, localGFEADOLCStiffnessScalar);
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Prepare the iterate x where we want to assemble - cos(identity) in 2D with z = x^2
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ auto deformationPowerBasis = makeBasis(
+ gridView,
+ power<gridDim>(
+ lagrange<2>()
+ ));
+ BlockVector<FieldVector<double,gridDim> > identity(compositeBasis.size({0}));
+ Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,gridDim> x){ return x; });
+ BlockVector<FieldVector<double,dim> > initialDeformation(compositeBasis.size({0}));
+ initialDeformation = 0;
+
+ Vector x;
+ RBMVector xRBM;
+ xRBM.resize(compositeBasis.size({0}));
+ x[_0].resize(compositeBasis.size({0}));
+ x[_1].resize(compositeBasis.size({1}));
+ for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < gridDim; j++)
+ initialDeformation[i][j] = std::cos(identity[i][j]);
+ initialDeformation[i][2] = identity[i][0]*identity[i][0];
+ x[_0][i] = initialDeformation[i];
+ xRBM[i].r = initialDeformation[i];
+ }
+
+ //////////////////////////////////////////////////////////////////////////////
+ // Assemble the Hessian using vector-mode and scalar-mode
+ //////////////////////////////////////////////////////////////////////////////
+
+ CorrectionTypeMixed gradientMixed;
+ gradientMixed[_0].resize(x[_0].size());
+ gradientMixed[_1].resize(x[_1].size());
+
+ MatrixTypeMixed hessianMatrixMixedScalar;
+ mixedAssemblerScalar.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixedScalar, true);
+
+ MatrixTypeMixed hessianMatrixMixedVector;
+ mixedAssemblerVector.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixedVector, true);
+
+ auto differenceMixed = hessianMatrixMixedScalar;
+ differenceMixed -= hessianMatrixMixedVector;
+
+ CorrectionType gradient;
+ gradient.resize(xRBM.size());
+
+ MatrixType hessianMatrixScalar;
+ assemblerScalar.assembleGradientAndHessian(xRBM, gradient, hessianMatrixScalar, true);
+
+ MatrixType hessianMatrixVector;
+ assemblerVector.assembleGradientAndHessian(xRBM, gradient, hessianMatrixVector, true);
+
+ auto difference = hessianMatrixScalar;
+ difference -= hessianMatrixVector;
+
+ if (differenceMixed.frobenius_norm() > 1e-8)
+ {
+ std::cerr << "MixedLocalGFEADOLCStiffness: The ADOL-C scalar mode and vector mode produce different Hessians!" << std::endl;
+ return 1;
+ }
+ if (difference.frobenius_norm() > 1e-8)
+ {
+ std::cerr << "LocalGFEADOLCStiffness: The ADOL-C scalar mode and vector mode produce different Hessians!" << std::endl;
+ return 1;
+ }
+ return 0;
+}