diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
index 50692f4f48ce7784f6cc67182c32bd52bbc8f4d9..c52d50af4bd69ce37b9169bed2fb3dc070a18f37 100644
--- a/dune/elasticity/assemblers/feassembler.hh
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -5,13 +5,15 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
-#include "localfestiffness.hh"
+// #include "localfestiffness.hh"
/** \brief A global FE assembler for variational problems
*/
-template <class Basis, class VectorType>
+template <class Basis, class VectorType, class LocalFEStiffness>
class FEAssembler {
typedef typename Basis::GridView GridView;
@@ -33,15 +35,13 @@ public:
protected:
- LocalFEStiffness<GridView,
- typename Basis::LocalFiniteElement,
- VectorType>* localStiffness_;
+ LocalFEStiffness* localStiffness_;
public:
/** \brief Constructor for a given grid */
FEAssembler(const Basis& basis,
- LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
+ LocalFEStiffness* localStiffness)
: basis_(basis),
localStiffness_(localStiffness)
{}
@@ -65,8 +65,8 @@ public:
-template <class Basis, class TargetSpace>
-void FEAssembler<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class LocalFEStiffness>
+void FEAssembler<Basis,TargetSpace,LocalFEStiffness>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
int n = basis_.size();
@@ -94,8 +94,8 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
}
-template <class Basis, class VectorType>
-void FEAssembler<Basis,VectorType>::
+template <class Basis, class VectorType, class LocalFEStiffness>
+void FEAssembler<Basis,VectorType,LocalFEStiffness>::
assembleGradientAndHessian(const VectorType& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
@@ -152,8 +152,8 @@ assembleGradientAndHessian(const VectorType& sol,
}
-template <class Basis, class VectorType>
-double FEAssembler<Basis, VectorType>::
+template <class Basis, class VectorType, class LocalFEStiffness>
+double FEAssembler<Basis, VectorType, LocalFEStiffness>::
computeEnergy(const VectorType& sol) const
{
double energy = 0;
@@ -179,4 +179,301 @@ computeEnergy(const VectorType& sol) const
return energy;
}
+
+
+
+
+
+
+
+/** \brief MultiTypeBlockvector extension of FEAssembler
+ * The Basis is assumed to be in correct shape
+ *
+ * Basis: We expect here a Dune::Functions basis
+ */
+template <class Basis, class BlockVectorType0, class BlockVectorType1, class LocalFEStiffness>
+class FEAssembler<Basis,Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>, LocalFEStiffness>
+{
+
+ using GridView = typename Basis::GridView;
+
+ //! Dimension of the grid.
+ enum { gridDim = GridView::dimension };
+
+ using VectorType = Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>;
+
+ using BT0 = typename BlockVectorType0::block_type;
+ using BT1 = typename BlockVectorType1::block_type;
+
+ using DT = typename BlockVectorType0::field_type;
+ static_assert(std::is_same<DT,typename BlockVectorType1::field_type>::value);
+
+ //!
+ using MatrixType = Dune::MultiTypeBlockMatrix<
+ Dune::MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<DT, BT0::dimension, BT0::dimension>>, Dune::BCRSMatrix<Dune::FieldMatrix<DT, BT0::dimension, BT1::dimension>> >,
+ Dune::MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<DT, BT1::dimension, BT0::dimension>>, Dune::BCRSMatrix<Dune::FieldMatrix<DT, BT1::dimension, BT1::dimension>> > >;
+
+public:
+ const Basis basis_;
+
+ //! Partitions to assemble the Gradient and Hessian, is interior here, because when doing parallel calculations, the assembled matrix and gradient are treated as if they do not overlap
+ static constexpr auto paritionsOverlap = Dune::Partitions::interiorBorder;
+
+protected:
+
+ std::shared_ptr<LocalFEStiffness> localStiffness_;
+
+public:
+
+ /** \brief Constructor for a given grid */
+ FEAssembler(const Basis& basis,
+ std::shared_ptr<LocalFEStiffness> localStiffness)
+ : basis_(basis),
+ localStiffness_(localStiffness)
+ {}
+
+ /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+ *
+ * This may be more efficient than computing them separately
+ */
+ virtual void assembleGradientAndHessian(const VectorType& sol,
+ VectorType& gradient,
+ MatrixType& hessian,
+ bool computeOccupationPattern=true) const;
+
+ /** \brief Compute the energy of a deformation state */
+ virtual double computeEnergy(const VectorType& sol) const;
+
+ //protected:
+ void getNeighborsPerVertex(std::array<std::array<Dune::MatrixIndexSet,2>,2>& nb) const;
+
+}; // end class
+
+
+
+template <class Basis, class BlockVectorType0, class BlockVectorType1, class LocalFEStiffness>
+void FEAssembler<Basis,Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>,LocalFEStiffness>::
+getNeighborsPerVertex(std::array<std::array<Dune::MatrixIndexSet,2>,2>& nb) const
+{
+ using namespace Dune::Indices;
+
+ // Set sizes of the 2x2 submatrices
+ for (size_t i=0; i<2; i++)
+ for (size_t j=0; j<2; j++)
+ {
+ nb[i][j].resize(basis_.size({i}), basis_.size({j}));
+ }
+
+ // A view on the FE basis on a single element
+ auto localView = basis_.localView();
+
+ // Loop over all leaf elements
+ for(const auto& e : elements(basis_.gridView()))
+ {
+ // Bind the local FE basis view to the current element
+ localView.bind(e);
+
+ // Add element stiffness matrix onto the global stiffness matrix
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ // The global index of the i-th local degree of freedom of the element 'e'
+ auto row = localView.index(i);
+
+ for (size_t j=0; j<localView.size(); j++ )
+ {
+ // The global index of the j-th local degree of freedom of the element 'e'
+ auto col = localView.index(j);
+
+ nb[row[0]][col[0]].add(row[1],col[1]);
+ }
+ }
+ }
+}
+
+template <class Basis, class BlockVectorType0, class BlockVectorType1, class LocalFEStiffness>
+void FEAssembler<Basis,Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>,LocalFEStiffness>::
+assembleGradientAndHessian(const VectorType& sol,
+ VectorType& gradient,
+ MatrixType& hessian, // type too long
+ bool computeOccupationPattern) const
+{
+ using namespace Dune::Indices;
+
+ if (computeOccupationPattern) {
+
+ std::array<std::array<Dune::MatrixIndexSet,2>,2> neighborsPerVertex;
+ getNeighborsPerVertex(neighborsPerVertex);
+ neighborsPerVertex[0][0].exportIdx(hessian[_0][_0]);
+ neighborsPerVertex[0][1].exportIdx(hessian[_0][_1]);
+ neighborsPerVertex[1][0].exportIdx(hessian[_1][_0]);
+ neighborsPerVertex[1][1].exportIdx(hessian[_1][_1]);
+
+ }
+
+ using namespace Dune::Indices;
+ using VectorType = Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>;
+
+ hessian = 0;
+
+ gradient[_0].resize(sol[_0].size());
+ gradient[_1].resize(sol[_1].size());
+ gradient = 0;
+
+ for (const auto& element : elements(basis_.gridView(), paritionsOverlap))
+ {
+ auto localView = basis_.localView();
+ localView.bind(element);
+
+ // size of all dof'S
+ auto flatSize0 = localView.tree().child(_0).size();
+ auto flatSize1 = localView.tree().child(_1).size();
+
+ // block sizes of the vector spaces
+ auto blockSize0 = BlockVectorType0::block_type::dimension;
+ auto blockSize1 = BlockVectorType1::block_type::dimension;
+
+ // entries in the two solution blocks
+ auto numOfBaseFct0 = flatSize0/blockSize0;
+ auto numOfBaseFct1 = flatSize1/blockSize1;
+
+ // Extract local solution
+ VectorType localSolution;
+ localSolution[_0].resize(numOfBaseFct0);
+ localSolution[_1].resize(numOfBaseFct1);
+
+ for (int i=0; i<numOfBaseFct0; i++)
+ {
+ // HACK
+ // get in global index of the i-th vertex of the element => this is the index of the P1 function
+ auto index = basis_.gridView().indexSet().subIndex(element,i,gridDim);
+ localSolution[_0][i] = sol[_0][index]; // vector copy!
+ }
+
+ if ( numOfBaseFct1 != 1 )
+ DUNE_THROW(Dune::Exception,"XXX");
+
+ // get in global index of the element => index of P0 elements
+ auto index = basis_.gridView().indexSet().index(element);
+ localSolution[_1][index] = sol[_1][index]; // vector copy!
+
+ VectorType localGradient;
+ localGradient[_0].resize(numOfBaseFct0);
+ localGradient[_1].resize(numOfBaseFct1);
+
+ // setup local matrix and gradient
+ const auto& lfe0 = localView.tree().child(_0,0).finiteElement();
+
+ localStiffness_->assembleGradientAndHessian(element, lfe0, localSolution, localGradient);
+
+ // Add element matrix to global stiffness matrix
+ for(int i=0; i<flatSize0; i++)
+ {
+ auto row = localView.index(i);
+
+ for (int j=0; j<flatSize0; j++ )
+ {
+ auto col = localView.index(j);
+ hessian[_0][_0][row[1]][col[1]][row[2]][col[2]] += localStiffness_->A_[i][j];
+ }
+ // continue counting
+ for (int j=flatSize0; j<flatSize0+flatSize1; j++ )
+ {
+ auto col = localView.index(j);
+ hessian[_0][_1][row[1]][col[1]][row[2]][col[2]] += localStiffness_->A_[i][j];
+ }
+ }
+ for(int i=flatSize0; i<flatSize0+flatSize1; i++)
+ {
+ auto row = localView.index(i);
+
+ for (int j=0; j<flatSize0; j++ )
+ {
+ auto col = localView.index(j);
+ hessian[_1][_0][row[1]][col[1]][row[2]][col[2]] += localStiffness_->A_[i][j];
+ }
+ for (int j=flatSize0; j<flatSize0+flatSize1; j++ )
+ {
+ auto col = localView.index(j);
+ hessian[_1][_1][row[1]][col[1]][row[2]][col[2]] += localStiffness_->A_[i][j];
+ }
+ }
+
+ // Add local gradient to global gradient
+ for (int i=0; i<numOfBaseFct0; i++)
+ {
+ // HACK
+ // get in global index of the i-th vertex of the element => this is the index of the P1 function
+ auto index = basis_.gridView().indexSet().subIndex(element,i,gridDim);
+ gradient[_0][index] += localGradient[_0][i];
+ }
+
+ if ( numOfBaseFct1 != 1 )
+ DUNE_THROW(Dune::Exception,"XXX");
+
+ // get in global index of the element => index of P0 elements
+ index = basis_.gridView().indexSet().index(element);
+ gradient[_1][index] += localGradient[_1][0];
+ }
+}
+
+
+template <class Basis, class BlockVectorType0, class BlockVectorType1, class LocalFEStiffness>
+double FEAssembler<Basis,Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>, LocalFEStiffness>::
+computeEnergy(const VectorType& sol) const
+{
+ using namespace Dune::Indices;
+ double energy = 0;
+
+ if (sol[_0].size()!=basis_.size({_0}) || sol[_1].size()!=basis_.size({_1}))
+ DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
+
+ auto localView = basis_.localView();
+
+ // Loop over all elements
+ for (const auto& element : elements(basis_.gridView(), paritionsOverlap))
+ {
+ localView.bind(element);
+
+ // size of all dof'S
+ auto flatSize0 = localView.tree().child(_0).size();
+ auto flatSize1 = localView.tree().child(_1).size();
+
+ // block sizes of the vector spaces
+ auto blockSize0 = BlockVectorType0::block_type::dimension;
+ auto blockSize1 = BlockVectorType1::block_type::dimension;
+
+ // entries in the two solution blocks
+ auto numOfBaseFct0 = flatSize0/blockSize0;
+ auto numOfBaseFct1 = flatSize1/blockSize1;
+
+ std::cerr << flatSize0 << " " << flatSize1 << " " << blockSize0 << " " << blockSize1 << std::endl;
+
+ // Extract local solution
+ VectorType localSolution;
+ localSolution[_0].resize(numOfBaseFct0);
+ localSolution[_1].resize(numOfBaseFct1);
+
+ for (int i=0; i<numOfBaseFct0; i++)
+ {
+ // get in global index of the i-th vertex of the element => this is the index of the P1 function
+ auto index = basis_.gridView().indexSet().subIndex(element,i,gridDim);
+ localSolution[_0][i] = sol[_0][index]; // vector copy!
+ }
+
+ if ( numOfBaseFct1 != 1 )
+ DUNE_THROW(Dune::Exception,"XXX");
+
+ // get in global index of the element => index of P0 elements
+ auto index = basis_.gridView().indexSet().index(element);
+ localSolution[_1][index] = sol[_1][index]; // vector copy!
+
+ const auto& lfe0 = localView.tree().child(_0,0).finiteElement();
+
+ energy += localStiffness_->energy(element, lfe0, localSolution);
+
+ }
+
+ return energy;
+
+}
#endif
diff --git a/dune/elasticity/assemblers/localadolcstiffness.hh b/dune/elasticity/assemblers/localadolcstiffness.hh
index 8afcfade68ff54ab8a8f659bd3c31c3accedb5c0..69c7f39261296438bea6f422ba4d3a5e23b278a8 100644
--- a/dune/elasticity/assemblers/localadolcstiffness.hh
+++ b/dune/elasticity/assemblers/localadolcstiffness.hh
@@ -10,10 +10,15 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/common/indices.hh>
+
#include <dune/istl/matrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/bvector.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
+
//#define ADOLC_VECTOR_MODE
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
@@ -195,4 +200,269 @@ assembleGradientAndHessian(const Entity& element,
#endif
}
+
+
+
+
+
+/** \brief Extension of LocalADOLCStiffness for MultiTypeVector formats
+ */
+template<class GV, class LFE, class BlockVectorType0, class BlockVectorType1>
+class LocalADOLCStiffness<GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>>
+ : public LocalFEStiffness<GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>>
+{
+ using GridView = GV;
+ using LocalFiniteElement = LFE;
+
+ // grid types
+ typedef typename GridView::Grid::ctype DT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+
+ // some other sizes
+ enum {gridDim=GridView::dimension};
+
+ using VectorType = Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>;
+
+ using BT0 = typename BlockVectorType0::block_type;
+ using BT1 = typename BlockVectorType1::block_type;
+
+ // HACK -> we convert the BlockTypes to same-looking types with value "abouble"
+ using AVectorType = Dune::MultiTypeBlockVector<Dune::BlockVector<Dune::FieldVector<adouble,BT0::dimension>>,
+ Dune::BlockVector<Dune::FieldVector<adouble,BT1::dimension>>>;
+
+public:
+
+ LocalADOLCStiffness(const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
+ : localEnergy_(energy)
+ {}
+
+ /** \brief Compute the energy at the current configuration */
+ virtual double energy (const Entity& e,
+ const LocalFiniteElement& localFiniteElement,
+ const VectorType& localConfiguration) const;
+
+ /** \brief Assemble the local stiffness matrix at the current position
+
+ This uses the automatic differentiation toolbox ADOL_C.
+ */
+ virtual void assembleGradientAndHessian(const Entity& e,
+ const LocalFiniteElement& localFiniteElement,
+ const VectorType& localConfiguration,
+ VectorType& localGradient);
+
+ const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
+
+};
+
+
+template<class GV, class LFE, class BlockVectorType0, class BlockVectorType1>
+double LocalADOLCStiffness<GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>>::
+energy(const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const VectorType& localSolution) const
+{
+ using namespace Dune::Indices;
+
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+ double pureEnergy;
+
+ using BT0 = typename BlockVectorType0::block_type;
+ using BT1 = typename BlockVectorType1::block_type;
+
+ // HACK -> we convert the BlockTypes to same-looking types with value "abouble"
+ using AVectorType = Dune::MultiTypeBlockVector<Dune::BlockVector<Dune::FieldVector<adouble,BT0::dimension>>,
+ Dune::BlockVector<Dune::FieldVector<adouble,BT1::dimension>>>;
+
+ AVectorType localASolution;
+ localASolution[_0].resize(localSolution[_0].size());
+ localASolution[_1].resize(localSolution[_1].size());
+
+
+ trace_on(rank);
+
+ adouble energy = 0;
+
+ for (size_t i=0; i<localSolution[_0].size(); i++)
+ for (size_t j=0; j<localSolution[_0][i].size(); j++)
+ localASolution[_0][i][j] <<= localSolution[_0][i][j];
+
+ for (size_t i=0; i<localSolution[_1].size(); i++)
+ for (size_t j=0; j<localSolution[_1][i].size(); j++)
+ localASolution[_1][i][j] <<= localSolution[_1][i][j];
+
+ energy = localEnergy_->energy(element,localFiniteElement,localASolution);
+
+ energy >>= pureEnergy;
+
+ trace_off(rank);
+
+ return pureEnergy;
+}
+
+
+
+// ///////////////////////////////////////////////////////////
+// Compute gradient and Hessian together
+// To compute the Hessian we need to compute the gradient anyway, so we may
+// as well return it. This saves assembly time.
+// ///////////////////////////////////////////////////////////
+template<class GV, class LFE, class BlockVectorType0, class BlockVectorType1>
+void LocalADOLCStiffness<GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>>::
+assembleGradientAndHessian(const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const VectorType& localSolution,
+ VectorType& localGradient)
+{
+ using namespace Dune::Indices;
+
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+ // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
+ energy(element, localFiniteElement, localSolution);
+
+ /////////////////////////////////////////////////////////////////
+ // Compute the energy gradient
+ /////////////////////////////////////////////////////////////////
+
+ using BT0 = typename BlockVectorType0::block_type;
+ using BT1 = typename BlockVectorType1::block_type;
+
+ // Compute the actual gradient
+ size_t nDofs0 = localSolution[_0].size();
+ size_t nDofs1 = localSolution[_1].size();
+ size_t blocksize0 = BT0::dimension;
+ size_t blocksize1 = BT1::dimension;
+
+ size_t nDoubles = nDofs0 * blocksize0 + nDofs1 * blocksize1;
+
+ std::vector<double> xp(nDoubles);
+ int idx=0;
+ for (size_t i=0; i<nDofs0; i++)
+ for (size_t j=0; j<blocksize0; j++)
+ xp[idx++] = localSolution[_0][i][j];
+
+ for (size_t i=0; i<nDofs1; i++)
+ for (size_t j=0; j<blocksize1; j++)
+ xp[idx++] = localSolution[_1][i][j];
+
+ // Compute gradient
+ std::vector<double> g(nDoubles);
+ gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
+
+ // Copy into Dune type
+ localGradient[_0].resize(localSolution[_0].size());
+ localGradient[_1].resize(localSolution[_1].size());
+
+ idx=0;
+ for (size_t i=0; i<nDofs0; i++)
+ for (size_t j=0; j<blocksize0; j++)
+ localGradient[_0][i][j] = g[idx++];
+
+ for (size_t i=0; i<nDofs1; i++)
+ for (size_t j=0; j<blocksize1; j++)
+ localGradient[_1][i][j] = g[idx++];
+
+ /////////////////////////////////////////////////////////////////
+ // Compute Hessian
+ /////////////////////////////////////////////////////////////////
+
+ // we use a flat matrix format for the local hessian
+ this->A_.setSize(nDoubles,nDoubles);
+
+#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<nDoubles; i++ )
+ {
+ // set the i-th unit vector
+ v[i] = 1.;
+
+ 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++)
+ this->A_[i][j] = w[j];
+
+ // reset the i-th unit vector
+ v[i] = 0.;
+ }
+
+
+// for (size_t i=0; i<nDofs0; i++)
+// for (size_t ii=0; ii<blocksize0; ii++)
+// {
+// // Evaluate Hessian in the direction of each vector of the orthonormal frame
+// for (size_t k=0; k<blocksize0; k++)
+// v[i*blocksize0 + k] = (k==ii);
+//
+// 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++)
+// this->A_[i*blocksize0 + ii][j] = w[j];
+//
+// // Make v the null vector again
+// std::fill(&v[i*blocksize0], &v[(i+1)*blocksize0], 0.0);
+// }
+// for (size_t i=0; i<nDofs1; i++)
+// for (size_t ii=0; ii<blocksize1; ii++)
+// {
+// // Evaluate Hessian in the direction of each vector of the orthonormal frame
+// for (size_t k=0; k<blocksize1; k++)
+// v[i*blocksize1 + k] = (k==ii);
+//
+// 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<nDofs0*blocksize0; j++)
+// this->A_[_1][_0][i][j/blocksize0][ii][j%blocksize0] = w[j];
+// for (size_t j=0; j<nDofs1*blocksize1; j++)
+// this->A_[_1][_1][i][j/blocksize1][ii][j%blocksize1] = w[j + nDofs0*blocksize0];
+//
+// // Make v the null vector again
+// std::fill(&v[i*blocksize1], &v[(i+1)*blocksize1], 0.0);
+// }
+#else
+ DUNE_THROW(Dune::NotImplemented,"this code should never be touched");
+// 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 (size_t j=0; j<nDirections; j++)
+// this->A_[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]);
+// }
+#endif
+}
+
#endif
diff --git a/dune/elasticity/assemblers/localenergy.hh b/dune/elasticity/assemblers/localenergy.hh
index 475ed7f340e276b67e4b0d2c50f3709a7badd54e..a5df849e110c14d6a1796024e5f4f5b2c3ab8704 100644
--- a/dune/elasticity/assemblers/localenergy.hh
+++ b/dune/elasticity/assemblers/localenergy.hh
@@ -11,11 +11,12 @@ namespace Elasticity {
template<class GridView, class LocalFiniteElement, class VectorType>
class LocalEnergy
{
- typedef typename VectorType::value_type::field_type RT;
typedef typename GridView::template Codim<0>::Entity Element;
public:
+ using RT = typename VectorType::value_type::field_type;
+
/** \brief Compute the energy
*
* \param element A grid element
@@ -29,6 +30,33 @@ public:
};
+
+/** \brief Extension of LocalEnergy for MultiTypeBlockVectors */
+template<class GV, class LFE, class BlockVectorType0, class BlockVectorType1>
+class LocalEnergy<GV, LFE, MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>>
+{
+ using Element = typename GV::template Codim<0>::Entity;
+ using VectorType = MultiTypeBlockVector<BlockVectorType0, BlockVectorType1>;
+
+ using RT = typename BlockVectorType0::value_type::field_type;
+ // check for same type
+ static_assert(std::is_same<RT,typename BlockVectorType1::value_type::field_type>::value);
+
+public:
+
+ /** \brief Compute the energy
+ * *
+ * \param element A grid element
+ * \param LocalFiniteElement A finite element on the reference element
+ * \param localConfiguration The coefficients of a FE function on the current element
+ */
+ /** \brief Compute the energy at the current configuration */
+ virtual RT energy (const Element& element,
+ const LFE& localFiniteElement,
+ const VectorType& localConfiguration) const = 0;
+
+};
+
} // namespace Elasticity
} // namespace Dune
diff --git a/dune/elasticity/assemblers/localfestiffness.hh b/dune/elasticity/assemblers/localfestiffness.hh
index f2d13cb54d0c3f5da2fdcfe30c70efa5bfa2ca17..82e568cb7f184727336913069a43f059f28691cd 100644
--- a/dune/elasticity/assemblers/localfestiffness.hh
+++ b/dune/elasticity/assemblers/localfestiffness.hh
@@ -3,6 +3,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
+#include <dune/istl/bvector.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
@@ -35,5 +36,46 @@ public:
};
+
+
+/** \brief Extension of LocalFEStiffness with MultiTypeBlockVector arguments
+ * It assembles into a Dune::MultiTypeBlockMatrix
+ * Currently, only available for exactly two blocks
+ * */
+template<class GV, class LFE, class BlockVectorType0, class BlockVectorType1>
+class LocalFEStiffness< GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0, BlockVectorType1> >
+: public Dune::Elasticity::LocalEnergy< GV, LFE, Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1> >
+{
+ // grid types
+ using Entity = typename GV::template Codim<0>::Entity;
+
+ using BT0 = typename BlockVectorType0::block_type;
+ using BT1 = typename BlockVectorType1::block_type;
+
+ using RT = typename BlockVectorType0::field_type;
+ static_assert(std::is_same<RT,typename BlockVectorType1::field_type>::value);
+
+ using VectorType = Dune::MultiTypeBlockVector<BlockVectorType0,BlockVectorType1>;
+
+ // we use a flat matrix type for best compatibility with Dune::Functions
+ using MatrixType = Dune::Matrix<Dune::FieldMatrix<RT,1,1> >;
+
+ // some other sizes
+ enum {gridDim=GV::dimension};
+
+public:
+
+ /** \brief Assemble the local gradient and tangent matrix at the current position
+ */
+ virtual void assembleGradientAndHessian(const Entity& e,
+ const LFE& localFiniteElement,
+ const VectorType& localConfiguration,
+ VectorType& localGradient) = 0;
+
+ // assembled data
+ MatrixType A_;
+
+};
+
#endif
diff --git a/dune/elasticity/materials/integralenergy.hh b/dune/elasticity/materials/integralenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..35e2e4c9db9dbbdbfd7ce314aedaa4b4deb0feaf
--- /dev/null
+++ b/dune/elasticity/materials/integralenergy.hh
@@ -0,0 +1,98 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_INTEGRALENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_INTEGRALENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+#include <dune/elasticity/assemblers/localenergy.hh>
+
+namespace Dune::Elasticity {
+
+template<class Basis, class EnergyDensity, class field_type=double>
+class IntegralEnergy
+ : public Elasticity::LocalEnergy<typename Basis::GridView,
+ typename Basis::LocalView::Tree::FiniteElement,
+ std::vector<Dune::FieldVector<field_type,Basis::GridView::dimension> > >
+{
+ // grid types
+ using GridView = typename Basis::GridView;
+ using DT = typename GridView::ctype;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+
+ // some other sizes
+ constexpr static auto gridDim = GridView::dimension;
+
+public:
+
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ */
+ IntegralEnergy(const EnergyDensity& energyDensity)
+ : energyDensity_(energyDensity)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ field_type energy (const Element& e,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
+
+ const EnergyDensity energyDensity_;
+};
+
+template <class Basis, class EnergyDensity, class field_type>
+field_type
+IntegralEnergy<Basis, EnergyDensity, field_type>::
+energy(const Element& element,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
+{
+ assert(element.type() == localFiniteElement.type());
+
+ field_type energy = 0;
+
+ // store gradients of shape functions and base functions
+ std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+ std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+ int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+ : localFiniteElement.localBasis().order() * gridDim;
+
+ const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+ for (const auto& qp : quad)
+ {
+ const DT integrationElement = element.geometry().integrationElement(qp.position());
+
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+
+ // Global position
+ auto x = element.geometry().global(qp.position());
+
+ // Get gradients of shape functions
+ localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+
+ // compute gradients of base functions
+ for (size_t i=0; i<gradients.size(); ++i)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ // Deformation gradient
+ FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
+ for (size_t i=0; i<gradients.size(); i++)
+ for (int j=0; j<gridDim; j++)
+ deformationGradient[j].axpy(localConfiguration[i][j], gradients[i]);
+
+ // Integrate the energy density
+ energy += qp.weight() * integrationElement * energyDensity_(x, deformationGradient);
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_INTEGRALENERGY_HH
+
diff --git a/dune/elasticity/materials/mooneyrivlinenergy.hh b/dune/elasticity/materials/mooneyrivlinenergy.hh
index bb40a824f207559d0ba81c5209ff2bc99939b99d..2cce10a9f550e956d7a285fb01dc5de908798e37 100644
--- a/dune/elasticity/materials/mooneyrivlinenergy.hh
+++ b/dune/elasticity/materials/mooneyrivlinenergy.hh
@@ -71,8 +71,80 @@ public:
mooneyrivlin_03,
mooneyrivlin_k;
std::string mooneyrivlin_energy;
+
+ field_type energyDensity(const Dune::FieldVector<field_type,gridDim>& x, const Dune::FieldMatrix<field_type, gridDim, gridDim>& derivative)
+ {
+ field_type strainEnergyCiarlet = 0;
+ field_type strainEnergy = 0;
+ field_type strainEnergyWithLog = 0;
+ field_type strainEnergyWithSquare = 0;
+
+ Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
+
+ for (int i=0; i<gridDim; i++) {
+ for (int j=0; j<gridDim; j++) {
+ for (int k=0; k<gridDim; k++)
+ C[i][j] += derivative[k][i] * derivative[k][j];
+ }
+ }
+
+ //////////////////////////////////////////////////////////
+ // Eigenvalues of FTF
+ //////////////////////////////////////////////////////////
+
+ // eigenvalues of C, i.e., squared singular values \sigma_i of F
+ Dune::FieldVector<field_type, dim> sigmaSquared;
+ FMatrixHelp::eigenValues(C, sigmaSquared);
+
+ field_type normFSquared = derivative.frobenius_norm2();
+ field_type detF = derivative.determinant();
+
+ field_type normFinvSquared = 0;
+
+ field_type c2Tilde = 0;
+ for (int i = 0; i < dim; i++) {
+ normFinvSquared += 1/sigmaSquared[i];
+ // compute D, which is the sum of the squared eigenvalues
+ for (int j = i+1; j < dim; j++)
+ c2Tilde += sigmaSquared[j]*sigmaSquared[i];
+ }
+
+ field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
+ // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
+ c2Tilde /= pow(detF, 4.0/dim);
+ field_type c2TildeMinus3 = c2Tilde - 3;
+ field_type detFMinus1 = detF - 1;
+
+ // Add the local energy density
+ strainEnergyCiarlet += (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF));
+ strainEnergy = 0;
+ strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
+ mooneyrivlin_01 * c2TildeMinus3 +
+ mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
+ mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 +
+ mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3 +
+ mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 +
+ mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 +
+ mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 +
+ mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3;
+ field_type logDetF = log(detF);
+ strainEnergyWithLog += ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
+ strainEnergyWithSquare += ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
+
+ std::cout << std::scientific;
+ if (mooneyrivlin_energy == "log") {
+ return strainEnergyWithLog;
+ } else if (mooneyrivlin_energy == "square") {
+ return strainEnergyWithSquare;
+ }
+ std::cout << std::fixed;
+
+ return strainEnergyCiarlet;
+ };
+
};
+
template <class GridView, class LocalFiniteElement, class field_type>
field_type
MooneyRivlinEnergy<GridView, LocalFiniteElement, field_type>::
@@ -80,6 +152,9 @@ energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
{
+ return 0;
+}
+#if 0
assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
@@ -160,7 +235,7 @@ energy(const Entity& element,
field_type detFMinus1 = detF - 1;
// Add the local energy density
- strainEnergyCiarlet += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
+ strainEnergyCiarlet += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF));
strainEnergy = 0;
strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
mooneyrivlin_01 * c2TildeMinus3 +
@@ -171,7 +246,7 @@ energy(const Entity& element,
mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 +
mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 +
mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3;
- field_type logDetF = std::log(detF);
+ field_type logDetF = log(detF);
strainEnergyWithLog += weight * ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
strainEnergyWithSquare += weight * ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
}
@@ -187,6 +262,27 @@ energy(const Entity& element,
return strainEnergyCiarlet;
}
+
+
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 657179223c07b3fdf2f6f29ffdb6cd4edec1530d..ca7377f22416ebccf161984b1b40679d55cef6bc 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -41,6 +41,7 @@
#include <dune/elasticity/materials/neohookeenergy.hh>
#include <dune/elasticity/materials/neumannenergy.hh>
#include <dune/elasticity/materials/sumenergy.hh>
+#include <dune/elasticity/materials/integralenergy.hh>
// grid dimension
const int dim = 3;