diff --git a/CHANGELOG.md b/CHANGELOG.md
index 93f9ddf1880a4ba9a5dd845abcabd872ce3b29f6..1c3b976cf2ba435ff3e88c4c95d595597b7f12a5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,9 @@
- Introduce class `LocalDensity` for material-specific implementations
- Introduce class `LocalIntegralEnergy` to work with the densities
+- Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element
## Deprecations and removals
- Redundant implementations of `LocalEnergy` classes are now replaced by combining `LocalDensity` and `LocalIntegralEnergy`
+- Local energies and FE assemblers with `dune-fufem` scalar bases are now deprecated
diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
index 9e0c82f59a912538dbc951844985d5d0b1c96ed6..62d1751d44ac6d68c35b59cec1e6567bb2f02bb3 100644
--- a/dune/elasticity/assemblers/feassembler.hh
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -1,19 +1,188 @@
#ifndef DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
#define DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
-#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/fufem/assemblers/istlbackend.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
+#include <dune/functions/functionspacebases/concepts.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+
#include "localfestiffness.hh"
-/** \brief A global FE assembler for variational problems
+namespace Dune::Elasticity {
+
+/** \brief A global FE assembler for variational problems for dune-functions global bases
*/
template <class Basis, class VectorType>
class FEAssembler {
+ using LocalView = typename Basis::LocalView;
+ using field_type = typename VectorType::field_type;
+
+ //! Dimension of the grid.
+ enum { gridDim = LocalView::GridView::dimension };
+
+public:
+
+ const Basis basis_;
+
+ /** \brief Partition type on which to assemble
+ *
+ * We want a fully nonoverlapping partition, and therefore need to skip ghost elements.
+ */
+ static constexpr auto partitionType = Dune::Partitions::interiorBorder;
+
+protected:
+
+ std::shared_ptr<LocalFEStiffness<LocalView,field_type>> localStiffness_;
+
+public:
+
+ /** \brief Constructor for a given basis and local assembler */
+ FEAssembler(const Basis& basis,
+ std::shared_ptr<LocalFEStiffness<LocalView, field_type>> 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
+ */
+ template<class MatrixType>
+ void assembleGradientAndHessian(const VectorType& configuration,
+ VectorType& gradient,
+ MatrixType& hessian,
+ bool computeOccupationPattern=true) const;
+
+ /** \brief Compute the energy of a deformation state */
+ auto computeEnergy(const VectorType& configuration) const;
+
+
+}; // end class
+
+
+template <class Basis, class VectorType>
+template <class MatrixType>
+void FEAssembler<Basis,VectorType>::
+assembleGradientAndHessian(const VectorType& configuration,
+ VectorType& gradient,
+ MatrixType& hessian,
+ bool computeOccupationPattern) const
+{
+ // create backends for multi index access
+ auto hessianBackend = Dune::Fufem::ISTLMatrixBackend(hessian);
+ auto configurationBackend = Dune::Functions::istlVectorBackend(configuration);
+ auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
+
+ if (computeOccupationPattern)
+ {
+ auto patternBuilder = hessianBackend.patternBuilder();
+ patternBuilder.resize(basis_,basis_);
+
+ auto localView = basis_.localView();
+
+ for (const auto& element : elements(basis_.gridView(), partitionType))
+ {
+ localView.bind(element);
+ for (size_t i=0; i<localView.size(); i++)
+ for (size_t j=0; j<localView.size(); j++)
+ patternBuilder.insertEntry( localView.index(i), localView.index(j) );
+ }
+ patternBuilder.setupMatrix();
+ }
+
+ gradientBackend.resize(basis_);
+
+ hessian = 0;
+ gradient = 0;
+
+ auto localView = basis_.localView();
+
+ for (const auto& element : elements(basis_.gridView(), partitionType))
+ {
+ localView.bind(element);
+ const auto size = localView.size();
+
+ // Extract local configuration
+ std::vector<field_type> localConfiguration(size);
+ std::vector<field_type> localGradient(size);
+ Matrix<field_type> localHessian;
+ localHessian.setSize(size, size);
+
+ for (int i=0; i<size; i++)
+ localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+ // setup local matrix and gradient
+ localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);
+
+ for(int i=0; i<size; i++)
+ {
+ auto row = localView.index(i);
+ gradientBackend[row] += localGradient[i];
+
+ for (int j=0; j<size; j++ )
+ {
+ auto col = localView.index(j);
+ // fufem ISTLBackend uses operator() for entry access
+ hessianBackend(row,col) += localHessian[i][j];
+ }
+ }
+ }
+}
+
+template <class Basis, class VectorType>
+auto FEAssembler<Basis,VectorType>::
+computeEnergy(const VectorType& configuration) const
+{
+ double energy = 0;
+
+ if (configuration.dim()!=basis_.dimension())
+ DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
+
+ auto localView = basis_.localView();
+ // fufem multi index access
+ auto configurationBackend = Functions::istlVectorBackend(configuration);
+
+ // Loop over all elements
+ for (const auto& element : elements(basis_.gridView(), partitionType))
+ {
+ localView.bind(element);
+ const auto size = localView.size();
+
+ std::vector<field_type> localConfiguration(size);
+
+ for (size_t i=0; i<size; i++)
+ localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+ energy += localStiffness_->energy(localView, localConfiguration);
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+
+
+namespace Dune
+{
+
+/** \brief A global FE assembler for variational problems (old fufem bases version)
+ */
+template <class Basis, class VectorType >
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.")
+FEAssembler
+{
+
typedef typename Basis::GridView GridView;
//! Dimension of the grid.
@@ -182,4 +351,7 @@ computeEnergy(const VectorType& sol) const
return energy;
}
+
+} // namespace Dune
+
#endif
diff --git a/dune/elasticity/assemblers/localadolcstiffness.hh b/dune/elasticity/assemblers/localadolcstiffness.hh
index a22bb83dcd930a288b5ef390f2879cbe4751c985..59fe3ea8769a277d73b731c2d4a99baa51e13f20 100644
--- a/dune/elasticity/assemblers/localadolcstiffness.hh
+++ b/dune/elasticity/assemblers/localadolcstiffness.hh
@@ -14,12 +14,130 @@
#include <dune/elasticity/assemblers/localfestiffness.hh>
-//#define ADOLC_VECTOR_MODE
+namespace Dune::Elasticity {
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
-template<class GridView, class LocalFiniteElement, class VectorType>
+template<class LocalView>
class LocalADOLCStiffness
+ : public LocalFEStiffness<LocalView,double>
+{
+ enum {gridDim=LocalView::GridView::dimension};
+
+public:
+
+ // accepts ADOL_C's "adouble" only
+ LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy)
+ : localEnergy_(energy)
+ {}
+
+ /** \brief Compute the energy at the current configuration */
+ virtual double energy (const LocalView& localView,
+ const std::vector<double>& localConfiguration) const;
+
+ /** \brief Assemble the local stiffness matrix at the current position
+ *
+ * This uses the automatic differentiation toolbox ADOL_C.
+ */
+ virtual void assembleGradientAndHessian(const LocalView& localView,
+ const std::vector<double>& localConfiguration,
+ std::vector<double>& localGradient,
+ Dune::Matrix<double>& localHessian);
+
+ std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> localEnergy_;
+
+};
+
+
+template<class LocalView>
+double LocalADOLCStiffness<LocalView>::
+energy(const LocalView& localView,
+ const std::vector<double>& localConfiguration) const
+{
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+ double pureEnergy;
+
+ std::vector<adouble> localAConfiguration(localConfiguration.size());
+
+ trace_on(rank);
+
+ adouble energy = 0;
+
+ for (size_t i=0; i<localConfiguration.size(); i++)
+ localAConfiguration[i] <<= localConfiguration[i];
+
+ energy = localEnergy_->energy(localView,localAConfiguration);
+
+ energy >>= pureEnergy;
+
+ trace_off(rank);
+
+ return pureEnergy;
+}
+
+
+template<class LocalView>
+void LocalADOLCStiffness<LocalView>::
+assembleGradientAndHessian(const LocalView& localView,
+ const std::vector<double>& localConfiguration,
+ std::vector<double>& localGradient,
+ Dune::Matrix<double>& localHessian
+ )
+{
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+ // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
+ energy(localView, localConfiguration);
+
+ /////////////////////////////////////////////////////////////////
+ // Compute the energy gradient
+ /////////////////////////////////////////////////////////////////
+
+ // Compute the actual gradient
+ size_t nDoubles = localConfiguration.size();
+
+ // Compute gradient
+ localGradient.resize(nDoubles);
+ gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+
+ /////////////////////////////////////////////////////////////////
+ // Compute Hessian
+ /////////////////////////////////////////////////////////////////
+
+ localHessian.setSize(nDoubles,nDoubles);
+
+ 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++)
+ {
+ // Evaluate Hessian multiplied with the i-th canonical basis vector (i.e.: the i-th column of the Hessian)
+ v[i] = 1.;
+
+ int rc= 3;
+ MINDEC(rc, hess_vec(rank, nDoubles, const_cast<double*>(localConfiguration.data()), v.data(), w.data())); // ADOL-C won't accept pointers with const values
+ if (rc < 0)
+ DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+ for (size_t j=0; j<nDoubles; j++)
+ localHessian[i][j] = w[j];
+
+ // Make v the null vector again
+ v[i] = 0.;
+ }
+
+ // TODO: implement ADOLC_VECTOR_MODE variant
+}
+
+} // namespace Dune::Elasticity
+
+
+/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
+ */
+template<class GridView, class LocalFiniteElement, class VectorType>
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!")
+LocalADOLCStiffness
: public LocalFEStiffness<GridView,LocalFiniteElement,VectorType>
{
// grid types
@@ -38,7 +156,7 @@ public:
//! Dimension of a tangent space
enum { blocksize = VectorType::value_type::dimension };
- LocalADOLCStiffness(const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
+ LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
: localEnergy_(energy)
{}
@@ -56,7 +174,7 @@ public:
const VectorType& localConfiguration,
VectorType& localGradient);
- const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
+ const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
};
diff --git a/dune/elasticity/assemblers/localenergy.hh b/dune/elasticity/assemblers/localenergy.hh
index 475ed7f340e276b67e4b0d2c50f3709a7badd54e..f4c56f3ecfa97bf26bce6522883796067c54948c 100644
--- a/dune/elasticity/assemblers/localenergy.hh
+++ b/dune/elasticity/assemblers/localenergy.hh
@@ -3,13 +3,33 @@
#include <vector>
-namespace Dune {
+namespace Dune::Elasticity {
-namespace Elasticity {
+/** \brief Base class for energies defined by integrating over one grid element */
+template<class LocalView, class field_type = double>
+class LocalEnergy
+{
+
+public:
+
+ /** \brief Compute the energy at the current configuration
+ *
+ * \param localView Container with current element and local function basis
+ * \param localConfiguration The coefficients of a FE function on the current element
+ */
+ virtual field_type energy (const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const = 0;
+
+};
+
+} // namespace Dune::Elasticity
+
+namespace Dune {
/** \brief Base class for energies defined by integrating over one grid element */
template<class GridView, class LocalFiniteElement, class VectorType>
-class LocalEnergy
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalEnergy with LocalView concept!")
+LocalEnergy
{
typedef typename VectorType::value_type::field_type RT;
typedef typename GridView::template Codim<0>::Entity Element;
@@ -29,8 +49,6 @@ public:
};
-} // namespace Elasticity
-
} // namespace Dune
#endif // DUNE_ELASTICITY_ASSEMBLERS_LOCALENERGY_HH
diff --git a/dune/elasticity/assemblers/localfestiffness.hh b/dune/elasticity/assemblers/localfestiffness.hh
index f2d13cb54d0c3f5da2fdcfe30c70efa5bfa2ca17..5aa1f982e88d3193cd9d9be0241f1d2fcff1dddd 100644
--- a/dune/elasticity/assemblers/localfestiffness.hh
+++ b/dune/elasticity/assemblers/localfestiffness.hh
@@ -6,9 +6,33 @@
#include <dune/elasticity/assemblers/localenergy.hh>
-template<class GridView, class LocalFiniteElement, class VectorType>
+namespace Dune::Elasticity {
+
+template<class LocalView, class RT = double>
class LocalFEStiffness
-: public Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, VectorType>
+: public Dune::Elasticity::LocalEnergy<LocalView, RT>
+{
+ enum {gridDim=LocalView::GridView::dimension};
+
+public:
+
+ /** \brief Assemble the local gradient and tangent matrix at the current position
+ */
+ virtual void assembleGradientAndHessian(const LocalView& localView,
+ const std::vector<RT>& localConfiguration,
+ std::vector<RT>& localGradient,
+ Dune::Matrix<RT>& localHessian
+ ) = 0;
+
+};
+
+} // namespace Dune::Elasticity
+
+
+template<class GridView, class LocalFiniteElement, class VectorType>
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalFEStiffness with LocalView concept!")
+LocalFEStiffness
+: public Dune::LocalEnergy<GridView, LocalFiniteElement, VectorType>
{
// grid types
typedef typename GridView::Grid::ctype DT;
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 432b73a029fe6abefed83cc548915a07a03d413f..5c7f54064a2d7bae8901a73f4f01262fccacb0fd 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -8,9 +8,11 @@
#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+#include <dune/fufem/assemblers/istlbackend.hh>
// Using a monotone multigrid as the inner solver
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
@@ -27,7 +29,10 @@
template <class BasisType, class VectorType>
void TrustRegionSolver<BasisType,VectorType>::
setup(const typename BasisType::GridView::Grid& grid,
- const FEAssembler<DuneFunctionsBasis<BasisType>, VectorType>* assembler,
+ std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+ Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
+ const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -65,7 +70,11 @@ setup(const typename BasisType::GridView::Grid& grid,
mgSetup_->ignore(std::const_pointer_cast<Dune::BitSetVector<blocksize> >(ignoreNodes_));
BasisType feBasis(grid.levelGridView(grid.maxLevel()));
- DuneFunctionsBasis<BasisType> basis(feBasis);
+
+ std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+ Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ BasisType,
+ DuneFunctionsBasis<BasisType> > basis(feBasis);
innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
innerMultigridStep_->mu_ = mu;
@@ -75,8 +84,11 @@ setup(const typename BasisType::GridView::Grid& grid,
if (grid.comm().rank()==0)
std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
#else
- DuneFunctionsBasis<BasisType> basis(grid.leafGridView());
+ std::conditional_t< // do we have a dune-functions basis?
+ Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ BasisType,
+ DuneFunctionsBasis<BasisType> > basis(grid.leafGridView());
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
@@ -117,15 +129,46 @@ setup(const typename BasisType::GridView::Grid& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
- OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> operatorAssembler(basis, basis);
- LaplaceAssembler<GridType,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> laplaceStiffness;
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+
+ using ScalarMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
ScalarMatrixType localA;
- operatorAssembler.assemble(laplaceStiffness, localA);
+ std::conditional_t< // do we have a dune-functions basis?
+ Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType>,
+ OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> > operatorAssembler(basis,basis);
+
+ if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+ {
+ using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+ // construct a Fufem Basis Assembler
+ auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement>();
+ // transform it to a dune-functions assembler
+ auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+ laplaceStiffness.assemble(element,
+ localMatrix,
+ trialLocalView.tree().child(0).finiteElement(),
+ ansatzLocalView.tree().child(0).finiteElement());
+ };
+
+
+ auto matrixBackend = Dune::Fufem::istlMatrixBackend(localA);
+ auto patternBuilder = matrixBackend.patternBuilder();
+
+ operatorAssembler.assembleBulkPattern(patternBuilder);
+ patternBuilder.setupMatrix();
+
+ operatorAssembler.assembleBulkEntries(matrixBackend, localAssembler);
+ }
+ else
+ {
+ LaplaceAssembler<GridType,
+ typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
+ typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> laplaceStiffness;
+
+ operatorAssembler.assemble(laplaceStiffness, localA);
+ }
ScalarMatrixType* A = new ScalarMatrixType(localA);
@@ -136,12 +179,37 @@ setup(const typename BasisType::GridView::Grid& grid,
// This will be used to monitor the gradient
// //////////////////////////////////////////////////////////////////////////////////////
- MassAssembler<GridType,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> massStiffness;
ScalarMatrixType localMassMatrix;
- operatorAssembler.assemble(massStiffness, localMassMatrix);
+ if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+ {
+ // construct a Fufem Basis Assembler
+ using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+ auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement>();
+ // transform it to a dune-functions assembler
+ auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+ massStiffness.assemble(element,
+ localMatrix,
+ trialLocalView.tree().child(0).finiteElement(),
+ ansatzLocalView.tree().child(0).finiteElement());
+ };
+
+ auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
+ auto massPatternBuilder = massMatrixBackend.patternBuilder();
+
+ operatorAssembler.assembleBulkPattern(massPatternBuilder);
+ massPatternBuilder.setupMatrix();
+
+ operatorAssembler.assembleBulkEntries(massMatrixBackend, localMassAssembler);
+ }
+ else
+ {
+ MassAssembler<GridType,
+ typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
+ typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> massStiffness;
+
+ operatorAssembler.assemble(massStiffness, localMassMatrix);
+ }
ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
@@ -151,10 +219,19 @@ setup(const typename BasisType::GridView::Grid& grid,
// ////////////////////////////////////////////////////////////
hessianMatrix_ = std::make_shared<MatrixType>();
- Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
- assembler_->getNeighborsPerVertex(indices);
- indices.exportIdx(*hessianMatrix_);
-
+ if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+ {
+ auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
+ auto hessianPatternBuilder = hessianBackend.patternBuilder();
+ operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
+ hessianPatternBuilder.setupMatrix();
+ }
+ else
+ {
+ Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+ assembler_->getNeighborsPerVertex(indices);
+ indices.exportIdx(*hessianMatrix_);
+ }
#if !HAVE_DUNE_PARMG
innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
@@ -174,9 +251,26 @@ setup(const typename BasisType::GridView::Grid& grid,
// On the lower grid levels a hierarchy of P1/Q1 spaces is used again.
////////////////////////////////////////////////////////////////////////
- typedef BasisType Basis;
- bool isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
+ using Basis = BasisType;
+ bool isP1Basis;
+ if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+ {
+ const auto dim = VectorType::value_type::dimension;
+ using namespace Dune::Functions::BasisFactory;
+ auto p1PowerBasis = makeBasis(
+ basis.gridView(),
+ power<dim>(
+ lagrange<1>()
+ ));
+
+ using P1PowerBasis = decltype(p1PowerBasis);
+ isP1Basis = std::is_same<Basis,P1PowerBasis>::value;
+ }
+ else
+ {
+ isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
+ }
using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
@@ -184,16 +278,31 @@ setup(const typename BasisType::GridView::Grid& grid,
// Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
if (not isP1Basis)
{
+ TransferOperatorType pkToP1TransferMatrix;
+
+ if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+ {
+ const auto dim = VectorType::value_type::dimension;
+ using namespace Dune::Functions::BasisFactory;
+ auto p1PowerBasis = makeBasis(
+ basis.gridView(),
+ power<dim>(
+ lagrange<1>()
+ ));
+ assembleGlobalBasisTransferMatrix(pkToP1TransferMatrix,p1PowerBasis,basis);
+ }
+ else
+ {
P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
-
- TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
- P1NodalBasis<typename GridType::LeafGridView,double>,
- DuneFunctionsBasis<BasisType> >(pkToP1TransferMatrix,p1Basis,basis);
- auto pkToP1 = std::make_shared< TruncatedCompressedMGTransfer<CorrectionType> >();
- transferOperators.back() = pkToP1;
- std::shared_ptr<TransferOperatorType> topTransferOperator = std::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
- pkToP1->setMatrix(topTransferOperator);
+ P1NodalBasis<typename GridType::LeafGridView,double>,
+ DuneFunctionsBasis<BasisType> >(pkToP1TransferMatrix,p1Basis,basis);
+ }
+
+ auto pkToP1 = std::make_shared< TruncatedCompressedMGTransfer<CorrectionType> >();
+ transferOperators.back() = pkToP1;
+ std::shared_ptr<TransferOperatorType> topTransferOperator = std::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+ pkToP1->setMatrix(topTransferOperator);
}
// Now the P1/Q1 restriction operators for the remaining levels
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index 2043fb3a0166481d32f6a615a55e12d223900b06..9856856bde42fc227c53abe687f8fd2b2a197add 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -3,8 +3,10 @@
#include <memory>
#include <vector>
+#include <type_traits>
#include <dune/common/bitsetvector.hh>
+#include <dune/common/concept.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
@@ -19,6 +21,8 @@
#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/functions/functionspacebases/concepts.hh>
+
#include <dune/elasticity/assemblers/feassembler.hh>
#if HAVE_DUNE_PARMG
@@ -35,13 +39,14 @@
#endif
-/** \brief Trust-region solver */
+/** \brief Trust-region solver */
template <class BasisType, class VectorType>
class TrustRegionSolver
: public IterativeSolver<VectorType,
Dune::BitSetVector<VectorType::value_type::dimension> >
{
typedef typename BasisType::GridView::Grid GridType;
+ using GridView = typename BasisType::GridView;
const static int blocksize = VectorType::value_type::dimension;
@@ -63,11 +68,21 @@ public:
TrustRegionSolver()
: IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(std::shared_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
- {}
+ {
+#if HAVE_DUNE_PARMG
+ // TODO currently, we cannot use PARMG together with dune-functions bases
+ static_assert( not Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() ,
+ "Currently, dune-parmg and dune-functions bases cannot be used together."
+ );
+#endif
+ }
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
- const FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* assembler,
+ std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+ Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
+ const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>*> assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -133,7 +148,10 @@ protected:
std::shared_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- const FEAssembler<DuneFunctionsBasis<BasisType>, VectorType>* assembler_;
+ std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+ Dune::models< Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
+ const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler_;
#if HAVE_DUNE_PARMG
/** \brief The solver for the inner problem */
diff --git a/dune/elasticity/materials/localintegralenergy.hh b/dune/elasticity/materials/localintegralenergy.hh
index 1ba6759d57845360700910e271e38a984bbbc32a..5b5df9b3a98a8feb18a898dafe54928ffdcd91b1 100644
--- a/dune/elasticity/materials/localintegralenergy.hh
+++ b/dune/elasticity/materials/localintegralenergy.hh
@@ -11,23 +11,107 @@
namespace Dune::Elasticity {
-template<class GridView, class LocalFiniteElement, class field_type=double>
+template<class LocalView, class field_type=double>
class LocalIntegralEnergy
- : public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension>>>
+ : public Elasticity::LocalEnergy<LocalView,field_type>
{
- // grid types
+ using GridView = typename LocalView::GridView;
using DT = typename GridView::Grid::ctype;
+
+ enum {gridDim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with a local energy density
+ */
+ LocalIntegralEnergy(const std::shared_ptr<LocalDensity<gridDim,field_type,DT>>& ld)
+ : localDensity_(ld)
+ {}
+
+ /** \brief Virtual destructor */
+ virtual ~LocalIntegralEnergy()
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ field_type energy (const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const;
+
+protected:
+ const std::shared_ptr<LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
+
+};
+
+template <class LocalView, class field_type>
+field_type
+LocalIntegralEnergy<LocalView, field_type>::
+energy(const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const
+{
+ // powerBasis: grab the finite element of the first child
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto& element = localView.element();
+
+ 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[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+ // Integrate the energy density
+ energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+
+// deprecated implementation in namespace Dune
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::LocalIntegralEnergy")
+LocalIntegralEnergy
+ : public LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension>>>
+{
using Entity = typename GridView::template Codim<0>::Entity;
+ using DT = typename GridView::Grid::ctype;
// some other sizes
enum {gridDim=GridView::dimension};
- enum {dim=GridView::dimension};
public:
/** \brief Constructor with a local energy density
*/
- LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type,DT>>& ld)
+ LocalIntegralEnergy(const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>>& ld)
: localDensity_(ld)
{}
@@ -41,7 +125,7 @@ public:
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
protected:
- const std::shared_ptr<LocalDensity<dim,field_type,DT>> localDensity_ = nullptr;
+ const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
};
@@ -49,10 +133,10 @@ protected:
template <class GridView, class LocalFiniteElement, class field_type>
field_type
-LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>::
+LocalIntegralEnergy<GridView, LocalFiniteElement, field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
+ const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
@@ -96,7 +180,7 @@ energy(const Entity& element,
return energy;
}
-} // namespace Dune::Elasticity
+} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
diff --git a/dune/elasticity/materials/neumannenergy.hh b/dune/elasticity/materials/neumannenergy.hh
index c1e11721a30b341677dc23f19396a41adec158c0..2bb416eec2262c1c203f5cbd69c88a8978019c5d 100644
--- a/dune/elasticity/materials/neumannenergy.hh
+++ b/dune/elasticity/materials/neumannenergy.hh
@@ -8,11 +8,99 @@
#include <dune/elasticity/assemblers/localenergy.hh>
+namespace Dune::Elasticity {
+
+template<class LocalView, class field_type=double>
+class NeumannEnergy
+ : public Elasticity::LocalEnergy<LocalView, field_type>
+{
+ using GridView = typename LocalView::GridView;
+ using DT = typename GridView::Grid::ctype;
+ using ctype = typename GridView::ctype;
+
+ enum {dim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ */
+ NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView>>& neumannBoundary,
+ const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,dim> >* neumannFunction)
+ : neumannBoundary_(neumannBoundary),
+ neumannFunction_(neumannFunction)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ field_type energy (const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const
+ {
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto& element = localView.element();
+
+ field_type energy = 0;
+
+ for (auto&& it : intersections(neumannBoundary_->gridView(), element)) {
+
+ if (not neumannBoundary_ or not neumannBoundary_->contains(it))
+ continue;
+
+ int quadOrder = localFiniteElement.localBasis().order();
+
+ const auto& quad = Dune::QuadratureRules<ctype, dim-1>::rule(it.type(), quadOrder);
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<ctype,dim>& quadPos = it.geometryInInside().global(quad[pt].position());
+
+ const auto integrationElement = it.geometry().integrationElement(quad[pt].position());
+
+ // The value of the local function
+ std::vector<Dune::FieldVector<ctype,1> > shapeFunctionValues;
+ localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+ Dune::FieldVector<field_type,dim> value(0);
+ for (size_t i=0; i<localFiniteElement.size(); i++)
+ for (int j=0; j<dim; j++)
+ value[j] += shapeFunctionValues[i] * localConfiguration[ localView.tree().child(j).localIndex(i) ];
+
+ // Value of the Neumann data at the current position
+ Dune::FieldVector<double,dim> neumannValue;
+
+ if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_))
+ dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
+ else
+ neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+
+ // Only translational dofs are affected by the Neumann force
+ for (size_t i=0; i<neumannValue.size(); i++)
+ energy += (neumannValue[i] * value[i]) * quad[pt].weight() * integrationElement;
+
+ }
+
+ }
+
+ return energy;
+ }
+
+private:
+ /** \brief The Neumann boundary */
+ const std::shared_ptr<BoundaryPatch<GridView>> neumannBoundary_;
+
+ /** \brief The function implementing the Neumann data */
+ const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,dim> >* neumannFunction_;
+};
+
+} // namespace Dune::Elasticity
+
+
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeumannEnergy
-: public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::NeumannEnergy")
+NeumannEnergy
+: public LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension> > >
{
// grid types
typedef typename GridView::ctype ctype;
@@ -92,6 +180,7 @@ private:
const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,dim> >* neumannFunction_;
};
-}
+} // namespace Dune
+
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
diff --git a/dune/elasticity/materials/sumenergy.hh b/dune/elasticity/materials/sumenergy.hh
index 685d51e6de6b1cb8afb9efc15910d8a953cdd835..209924678da7438f649f412c95ab7befa7fc2469 100644
--- a/dune/elasticity/materials/sumenergy.hh
+++ b/dune/elasticity/materials/sumenergy.hh
@@ -11,16 +11,56 @@
#include <dune/elasticity/assemblers/localenergy.hh>
+namespace Dune::Elasticity {
+
+template<class LocalView, class field_type=double>
+class SumEnergy
+ : public Elasticity::LocalEnergy<LocalView, field_type>
+{
+ using GridView = typename LocalView::GridView;
+ using DT = typename GridView::Grid::ctype;
+
+ enum {dim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ */
+ SumEnergy(std::shared_ptr<Elasticity::LocalEnergy<LocalView, field_type>> a,
+ std::shared_ptr<Elasticity::LocalEnergy<LocalView, field_type>> b)
+ : a_(a),
+ b_(b)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ field_type energy (const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const
+ {
+ return a_->energy(localView, localConfiguration)
+ + b_->energy(localView, localConfiguration);
+ }
+
+private:
+ std::shared_ptr<Elasticity::LocalEnergy<LocalView, field_type>> a_;
+ std::shared_ptr<Elasticity::LocalEnergy<LocalView, field_type>> b_;
+};
+
+} // namespace Dune::Elasticity
+
+
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class SumEnergy
-: public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::SumEnergy")
+SumEnergy
+: public LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
{
// grid types
typedef typename GridView::ctype ctype;
typedef typename GridView::template Codim<0>::Entity Entity;
+ // some other sizes
enum {dim=GridView::dimension};
public:
@@ -28,8 +68,8 @@ public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
- SumEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > a,
- std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > b)
+ SumEnergy(std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > a,
+ std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > b)
: a_(a),
b_(b)
{}
@@ -45,11 +85,13 @@ public:
private:
- std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > a_;
+ std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > a_;
+
+ std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > b_;
- std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > b_;
};
-}
+} // namespace Dune
+
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_SUMENERGY_HH
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 848c1e5624fdbc5f5a0fa6a4e820bb92efb7541a..adf18afa768097fd76f4fb08e8b0ba015943e430 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -25,7 +25,6 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -151,16 +150,28 @@ int main (int argc, char *argv[]) try
GridView gridView = grid->leafGridView();
#endif
- // FE basis spanning the FE space that we are working in
- typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
- FEBasis feBasis(gridView);
+
+ ////////////////////////////////////////////
+ // Basis
+ ////////////////////////////////////////////
+
+ using namespace Dune::Functions::BasisFactory;
+ auto basis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ using PowerBasis = decltype(basis);
+ using LocalView = typename PowerBasis::LocalView;
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
+ BitSetVector<dim> dirichletVertices(gridView.size(dim), false);
+ BitSetVector<dim> neumannVertices(gridView.size(dim), false);
const GridView::IndexSet& indexSet = gridView.indexSet();
@@ -185,19 +196,19 @@ int main (int argc, char *argv[]) try
}
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
- std::cout << "On process " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+ std::cout << "On process " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
- BitSetVector<1> dirichletNodes(feBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+ BitSetVector<dim> dirichletNodes(basis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,basis,dirichletNodes);
- BitSetVector<1> neumannNodes(feBasis.size(), false);
- constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+ BitSetVector<dim> neumannNodes(basis.size(), false);
+ constructBoundaryDofs(*neumannBoundary,basis,neumannNodes);
- BitSetVector<dim> dirichletDofs(feBasis.size(), false);
- for (size_t i=0; i<feBasis.size(); i++)
+ BitSetVector<dim> dirichletDofs(basis.size(), false);
+ for (size_t i=0; i<basis.size(); i++)
if (dirichletNodes[i][0])
for (int j=0; j<dim; j++)
dirichletDofs[i][j] = true;
@@ -206,20 +217,12 @@ int main (int argc, char *argv[]) try
// Initial iterate
// //////////////////////////
- using namespace Dune::Functions::BasisFactory;
- auto powerBasis = makeBasis(
- gridView,
- power<dim>(
- lagrange<order>(),
- blockedInterleaved()
- ));
-
- SolutionType x(feBasis.size());
+ SolutionType x(basis.size());
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
- Dune::Functions::interpolate(powerBasis, x, pythonInitialDeformation);
+ Dune::Functions::interpolate(basis, x, pythonInitialDeformation);
////////////////////////////////////////////////////////
// Main homotopy loop
@@ -227,12 +230,12 @@ int main (int argc, char *argv[]) try
// Output initial iterate (of homotopy loop)
SolutionType identity;
- Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+ Dune::Functions::interpolate(basis, identity, [](FieldVector<double,dim> x){ return x; });
SolutionType displacement = x;
displacement -= identity;
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis, displacement);
auto localDisplacementFunction = localFunction(displacementFunction);
// We need to subsample, because VTK cannot natively display real second-order functions
@@ -288,33 +291,21 @@ int main (int argc, char *argv[]) try
if(!elasticDensity)
DUNE_THROW(Exception, "Error: Selected energy not available!");
- auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble>>(elasticDensity);
-
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(&neumannBoundary,neumannFunction.get());
+ auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
- SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> totalEnergy(elasticEnergy, neumannEnergy);
+ auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,adouble>>(neumannBoundary,neumannFunction.get());
- LocalADOLCStiffness<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- SolutionType> localADOLCStiffness(&totalEnergy);
+ auto totalEnergy = std::make_shared<Elasticity::SumEnergy<LocalView,adouble>>(elasticEnergy, neumannEnergy);
- // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFEBasis(feBasis);
+ auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(totalEnergy);
- FEAssembler<FufemFEBasis,SolutionType> assembler(fufemFEBasis, &localADOLCStiffness);
+ Elasticity::FEAssembler<PowerBasis,SolutionType> assembler(basis, localADOLCStiffness);
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
- TrustRegionSolver<FEBasis,SolutionType> solver;
+ TrustRegionSolver<PowerBasis,SolutionType> solver;
solver.setup(*grid,
&assembler,
x,
@@ -346,7 +337,7 @@ int main (int argc, char *argv[]) try
// Extract object member functions as Dune functions
PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
- Dune::Functions::interpolate(powerBasis, x, dirichletValues, dirichletDofs);
+ Dune::Functions::interpolate(basis, x, dirichletValues, dirichletDofs);
// /////////////////////////////////////////////////////
// Solve!
@@ -365,7 +356,7 @@ int main (int argc, char *argv[]) try
auto displacement = x;
displacement -= identity;
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));