diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 4e5fdbc0ecc069ef52a5b62d7e7b03580f621f91..8f7c34528b63bd2aa66b62c272ec8c5dd58acc24 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -37,9 +37,9 @@ class LocalFiniteElementFactory
public:
static auto get(const typename Basis::LocalView& localView,
std::integral_constant<std::size_t, i> iType)
- -> decltype(localView.tree().child(iType).finiteElement())
+ -> decltype(localView.tree().child(iType,0).finiteElement())
{
- return localView.tree().child(iType).finiteElement();
+ return localView.tree().child(iType,0).finiteElement();
}
};
diff --git a/dune/gfe/mixedgfeassembler.hh b/dune/gfe/mixedgfeassembler.hh
index 941b035f465d76b29e9fb030997a7c45b122dae2..2cd87af561b0a65813941a8ced4779eedfde5ab3 100644
--- a/dune/gfe/mixedgfeassembler.hh
+++ b/dune/gfe/mixedgfeassembler.hh
@@ -31,11 +31,9 @@ class MixedGFEAssembler {
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize0> > MatrixBlock10;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > MatrixBlock11;
+public:
typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixBlock00,MatrixBlock01>,
Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > MatrixType;
-
-protected:
-public:
const Basis basis_;
MixedLocalGeodesicFEStiffness<Basis,
@@ -195,26 +193,38 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
#endif
using namespace Dune::TypeTree::Indices;
- const int nDofs0 = localView.tree().child(_0).finiteElement().size();
- const int nDofs1 = localView.tree().child(_1).finiteElement().size();
+ const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
+ const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
+ // This loop reads out the pattern for a local matrix; in each element, we have localView.size() degrees of freedom; from the composite and powerbasis layers
+ // nDofs0 are the degrees of freedom for *one* subspacebasis of the power basis of the displacement part;
+ // nDofs1 are the degrees of freedom for *one* subspacebasis of the power basis of the rotational part
+ // this is why the indices (_0,0) and (_1,0) are used: _0 takes the whole displacement part and _1 the whole rotational part; and 0 the first subspacebasis respectively
// Extract local solution
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (int i=0; i<nDofs0+nDofs1; i++)
{
+ int localIndexI = 0;
+ if (i < nDofs0) {
+ auto& node = localView.tree().child(_0,0);
+ localIndexI = node.localIndex(i);
+ } else {
+ auto& node = localView.tree().child(_1,0);
+ localIndexI = node.localIndex(i-nDofs0);
+ }
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- if (localIndexSet.index(i)[0] == 0)
- localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
- else
- localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
+ auto multiIndex = localIndexSet.index(localIndexI);
#else
- if (localView.index(i)[0] == 0)
- localConfiguration0[i] = configuration0[localView.index(i)[1]];
- else
- localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
-#endif
+ auto multiIndex = localView.index(localIndexI);
+#endif
+ //CompositeBasis number is contained in multiIndex[0], the Subspacebasis is contained in multiIndex[2]
+ //multiIndex[1] contains the actual index
+ if (multiIndex[0] == 0)
+ localConfiguration0[i] = configuration0[multiIndex[1]];
+ else if (multiIndex[0] == 1)
+ localConfiguration1[i-nDofs0] = configuration1[multiIndex[1]];
}
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
@@ -228,18 +238,34 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0+nDofs1; i++)
{
+ int localIndexRow = 0;
+ if (i < nDofs0) {
+ auto& node = localView.tree().child(_0,0);
+ localIndexRow = node.localIndex(i);
+ } else {
+ auto& node = localView.tree().child(_1,0);
+ localIndexRow = node.localIndex(i-nDofs0);
+ }
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- auto row = localIndexSet.index(i);
+ auto row = localIndexSet.index(localIndexRow);
#else
- auto row = localView.index(i);
+ auto row = localView.index(localIndexRow);
#endif
for (int j=0; j<nDofs0+nDofs1; j++ )
{
+ int localIndexCol = 0;
+ if (j < nDofs0) {
+ auto& node = localView.tree().child(_0,0);
+ localIndexCol = node.localIndex(j);
+ } else {
+ auto& node = localView.tree().child(_1,0);
+ localIndexCol = node.localIndex(j-nDofs0);
+ }
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- auto col = localIndexSet.index(j);
+ auto col = localIndexSet.index(localIndexCol);
#else
- auto col = localView.index(j);
+ auto col = localView.index(localIndexCol);
#endif
if (row[0]==0 and col[0]==0)
@@ -256,17 +282,10 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
}
// Add local gradient to global gradient
-#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- if (localIndexSet.index(i)[0] == 0)
- gradient0[localIndexSet.index(i)[1]] += localGradient0[i];
- else
- gradient1[localIndexSet.index(i)[1]] += localGradient1[i-nDofs0];
-#else
- if (localView.index(i)[0] == 0)
- gradient0[localView.index(i)[1]] += localGradient0[i];
+ if (row[0] == 0)
+ gradient0[row[1]] += localGradient0[i];
else
- gradient1[localView.index(i)[1]] += localGradient1[i-nDofs0];
-#endif
+ gradient1[row[1]] += localGradient1[i-nDofs0];
}
}
@@ -343,25 +362,33 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
// Number of degrees of freedom on this element
using namespace Dune::TypeTree::Indices;
- const int nDofs0 = localView.tree().child(_0).finiteElement().size();
- const int nDofs1 = localView.tree().child(_1).finiteElement().size();
+ const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
+ const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (int i=0; i<nDofs0+nDofs1; i++)
{
+ int localIndexI = 0;
+ if (i < nDofs0) {
+ auto& node = localView.tree().child(_0,0);
+ localIndexI = node.localIndex(i);
+ } else {
+ auto& node = localView.tree().child(_1,0);
+ localIndexI = node.localIndex(i-nDofs0);
+ }
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- if (localIndexSet.index(i)[0] == 0)
- localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
- else
- localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
+ auto multiIndex = localIndexSet.index(localIndexI);
#else
- if (localView.index(i)[0] == 0)
- localConfiguration0[i] = configuration0[localView.index(i)[1]];
- else
- localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
-#endif
+ auto multiIndex = localView.index(localIndexI);
+#endif
+ // The CompositeBasis number is contained in multiIndex[0]
+ // multiIndex[1] contains the actual index
+ if (multiIndex[0] == 0)
+ localConfiguration0[i] = configuration0[multiIndex[1]];
+ else if (multiIndex[0] == 1)
+ localConfiguration1[i-nDofs0] = configuration1[multiIndex[1]];
}
energy += localStiffness_->energy(localView,
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 2cd74a7ee6a9956e54588af84fdc364137a4cad5..6a88c89460cd54351c55772c4505faa7b7dcfd65 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -487,7 +487,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
try {
energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
} catch (Dune::Exception &e) {
- std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+ std::cerr << "Error while computing the energy of the new iterate: " << e << std::endl;
std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
newIterate = x_;
solvedByInnerSolver = false;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 1bd4144a3c9a616571fffa131ce13bb40723b2ed..6d1719cf94d6fd103862869df0eb20bebbbea207 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -218,13 +218,17 @@ int main (int argc, char *argv[]) try
using namespace Dune::Functions::BasisFactory;
#ifdef MIXED_SPACE
+ const int dimRotation = Rotation<double,dim>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
composite(
- lagrange<displacementOrder>(),
- lagrange<rotationOrder>()
- )
- );
+ power<dim>(
+ lagrange<displacementOrder>()
+ ),
+ power<dimRotation>(
+ lagrange<rotationOrder>()
+ )
+ ));
typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;