diff --git a/dune/gfe/mixedcosseratenergy.hh b/dune/gfe/mixedcosseratenergy.hh
index 4239ede37e3d9beb851838e5dcb83dd259fab7a9..8c0bcf455dc7d4a4f31f402adab7ed1361087bc8 100644
--- a/dune/gfe/mixedcosseratenergy.hh
+++ b/dune/gfe/mixedcosseratenergy.hh
@@ -20,14 +20,16 @@
//#define QUADRATIC_MEMBRANE_ENERGY
-template<class GridView, class DisplacementLocalFiniteElement, class OrientationLocalFiniteElement, int dim, class field_type=double>
+template<class DisplacementBasis, class OrientationBasis, int dim, class field_type=double>
class MixedCosseratEnergy
- : public MixedLocalGeodesicFEStiffness<GridView,
- DisplacementLocalFiniteElement,RealTuple<field_type,dim>,
- OrientationLocalFiniteElement,Rotation<field_type,dim> >
+ : public MixedLocalGeodesicFEStiffness<DisplacementBasis,RealTuple<field_type,dim>,
+ OrientationBasis,Rotation<field_type,dim> >
{
// grid types
- typedef typename GridView::Grid::ctype DT;
+ typedef typename DisplacementBasis::LocalView::Tree::FiniteElement DisplacementLocalFiniteElement;
+ typedef typename OrientationBasis::LocalView::Tree::FiniteElement OrientationLocalFiniteElement;
+ typedef typename DisplacementBasis::GridView GridView;
+ typedef typename GridView::ctype DT;
typedef field_type RT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -263,11 +265,11 @@ public:
const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction_;
};
-template <class GridView, class DeformationLocalFiniteElement, class OrientationLocalFiniteElement, int dim, class field_type>
-typename MixedCosseratEnergy<GridView,DeformationLocalFiniteElement,OrientationLocalFiniteElement,dim,field_type>::RT
-MixedCosseratEnergy<GridView,DeformationLocalFiniteElement,OrientationLocalFiniteElement,dim,field_type>::
+template <class DeformationBasis, class OrientationBasis, int dim, class field_type>
+typename MixedCosseratEnergy<DeformationBasis,OrientationBasis,dim,field_type>::RT
+MixedCosseratEnergy<DeformationBasis,OrientationBasis,dim,field_type>::
energy(const Entity& element,
- const DeformationLocalFiniteElement& deformationLocalFiniteElement,
+ const DisplacementLocalFiniteElement& deformationLocalFiniteElement,
const std::vector<RealTuple<field_type,dim> >& localDeformationConfiguration,
const OrientationLocalFiniteElement& orientationLocalFiniteElement,
const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const
@@ -278,7 +280,7 @@ energy(const Entity& element,
RT energy = 0;
- typedef LocalGeodesicFEFunction<gridDim, DT, DeformationLocalFiniteElement, RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
+ typedef LocalGeodesicFEFunction<gridDim, DT, DisplacementLocalFiniteElement, RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
typedef LocalGeodesicFEFunction<gridDim, DT, OrientationLocalFiniteElement, Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
@@ -363,20 +365,20 @@ energy(const Entity& element,
if (not neumannFunction_)
return energy;
- for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
-
- if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
+ for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+ {
+ if (not neumannBoundary_ or not neumannBoundary_->contains(it))
continue;
const Dune::QuadratureRule<DT, gridDim-1>& quad
- = Dune::QuadratureRules<DT, gridDim-1>::rule(it->type(), quadOrder);
+ = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
+ const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
- const DT integrationElement = it->geometry().integrationElement(quad[pt].position());
+ const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function
RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
@@ -387,7 +389,7 @@ energy(const Entity& element,
if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
else
- neumannFunction_->evaluate(it->geometry().global(quad[pt].position()), neumannValue);
+ 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++)
@@ -400,5 +402,5 @@ energy(const Entity& element,
return energy;
}
-#endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
+#endif //#ifndef DUNE_GFE_MIXEDCOSSERATENERGY_HH
diff --git a/dune/gfe/mixedgfeassembler.hh b/dune/gfe/mixedgfeassembler.hh
index cfa4bc65a44dc8988d770c8ea64199797fbee1d6..ae917ea165044c01d8a33f15d8950ae538954a9c 100644
--- a/dune/gfe/mixedgfeassembler.hh
+++ b/dune/gfe/mixedgfeassembler.hh
@@ -35,20 +35,16 @@ public:
const Basis0 basis0_;
const Basis1 basis1_;
- MixedLocalGeodesicFEStiffness<GridView,
- typename Basis0::LocalFiniteElement,
- TargetSpace0,
- typename Basis1::LocalFiniteElement,
- TargetSpace1>* localStiffness_;
+ MixedLocalGeodesicFEStiffness<Basis0, TargetSpace0,
+ Basis1, TargetSpace1>* localStiffness_;
public:
/** \brief Constructor for a given grid */
MixedGFEAssembler(const Basis0& basis0,
const Basis1& basis1,
- MixedLocalGeodesicFEStiffness<GridView,
- typename Basis0::LocalFiniteElement, TargetSpace0,
- typename Basis0::LocalFiniteElement, TargetSpace1>* localStiffness)
+ MixedLocalGeodesicFEStiffness<Basis0, TargetSpace0,
+ Basis1, TargetSpace1>* localStiffness)
: basis0_(basis0),
basis1_(basis1),
localStiffness_(localStiffness)
@@ -94,50 +90,57 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
Dune::MatrixIndexSet& nb10,
Dune::MatrixIndexSet& nb11) const
{
- nb00.resize(basis0_.size(), basis0_.size());
- nb01.resize(basis0_.size(), basis1_.size());
- nb10.resize(basis1_.size(), basis0_.size());
- nb11.resize(basis1_.size(), basis1_.size());
+ nb00.resize(basis0_.indexSet().size(), basis0_.indexSet().size());
+ nb01.resize(basis0_.indexSet().size(), basis1_.indexSet().size());
+ nb10.resize(basis1_.indexSet().size(), basis0_.indexSet().size());
+ nb11.resize(basis1_.indexSet().size(), basis1_.indexSet().size());
+
+ // A view on the FE basis on a single element
+ typename Basis0::LocalView localView0(&basis0_);
+ typename Basis1::LocalView localView1(&basis1_);
+ auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+ auto localIndexSet1 = basis1_.indexSet().localIndexSet();
// Grid view must be the same for both bases
- ElementIterator it = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endit = basis0_.getGridView().template end<0,Dune::Interior_Partition> ();
+ ElementIterator it = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+ ElementIterator endit = basis0_.gridView().template end<0,Dune::Interior_Partition> ();
for (; it!=endit; ++it) {
- const typename Basis0::LocalFiniteElement& lfe0 = basis0_.getLocalFiniteElement(*it);
- const typename Basis1::LocalFiniteElement& lfe1 = basis1_.getLocalFiniteElement(*it);
+ // Bind the local FE basis view to the current element
+ localView0.bind(*it);
+ localView1.bind(*it);
+ localIndexSet0.bind(localView0);
+ localIndexSet1.bind(localView1);
- for (size_t i=0; i<lfe0.localBasis().size(); i++) {
+ for (size_t i=0; i<localView0.size(); i++) {
- int iIdx = basis0_.index(*it,i);
+ int iIdx = localIndexSet0.index(i)[0];
- for (size_t j=0; j<lfe0.localBasis().size(); j++) {
- int jIdx = basis0_.index(*it,j);
+ for (size_t j=0; j<localView0.size(); j++) {
+ int jIdx = localIndexSet0.index(j)[0];
nb00.add(iIdx, jIdx);
}
- for (size_t j=0; j<lfe1.localBasis().size(); j++) {
- int jIdx = basis1_.index(*it,j);
+ for (size_t j=0; j<localView1.size(); j++) {
+ int jIdx = localIndexSet1.index(j)[0];
nb01.add(iIdx, jIdx);
-
}
}
- for (size_t i=0; i<lfe1.localBasis().size(); i++) {
+ for (size_t i=0; i<localView1.size(); i++) {
- int iIdx = basis1_.index(*it,i);
+ int iIdx = localIndexSet1.index(i)[0];
- for (size_t j=0; j<lfe0.localBasis().size(); j++) {
- int jIdx = basis0_.index(*it,j);
+ for (size_t j=0; j<localView0.size(); j++) {
+ int jIdx = localIndexSet0.index(j)[0];
nb10.add(iIdx, jIdx);
}
- for (size_t j=0; j<lfe1.localBasis().size(); j++) {
- int jIdx = basis1_.index(*it,j);
+ for (size_t j=0; j<localView1.size(); j++) {
+ int jIdx = localIndexSet1.index(j)[0];
nb11.add(iIdx, jIdx);
-
}
}
@@ -184,70 +187,82 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
gradient1.resize(configuration1.size());
gradient1 = 0;
- ElementIterator it = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endit = basis0_.getGridView().template end<0,Dune::Interior_Partition> ();
+ // A view on the FE basis on a single element
+ typename Basis0::LocalView localView0(&basis0_);
+ typename Basis1::LocalView localView1(&basis1_);
+ auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+ auto localIndexSet1 = basis1_.indexSet().localIndexSet();
+
+ ElementIterator it = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+ ElementIterator endit = basis0_.gridView().template end<0,Dune::Interior_Partition> ();
for( ; it != endit; ++it ) {
- const int nDofs0 = basis0_.getLocalFiniteElement(*it).localBasis().size();
- const int nDofs1 = basis1_.getLocalFiniteElement(*it).localBasis().size();
+ // Bind the local FE basis view to the current element
+ localView0.bind(*it);
+ localView1.bind(*it);
+ localIndexSet0.bind(localView0);
+ localIndexSet1.bind(localView1);
+
+ const int nDofs0 = localView0.size();
+ const int nDofs1 = localView1.size();
// Extract local solution
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (int i=0; i<nDofs0; i++)
- localConfiguration0[i] = configuration0[basis0_.index(*it,i)];
+ localConfiguration0[i] = configuration0[localIndexSet0.index(i)[0]];
for (int i=0; i<nDofs1; i++)
- localConfiguration1[i] = configuration1[basis1_.index(*it,i)];
+ localConfiguration1[i] = configuration1[localIndexSet1.index(i)[0]];
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
std::vector<Dune::FieldVector<double,blocksize1> > localGradient1(nDofs1);
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(*it,
- basis0_.getLocalFiniteElement(*it), localConfiguration0,
- basis1_.getLocalFiniteElement(*it), localConfiguration1,
+ localView0.tree().finiteElement(), localConfiguration0,
+ localView1.tree().finiteElement(), localConfiguration1,
localGradient0, localGradient1);
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0; i++) {
- int row = basis0_.index(*it,i);
+ int row = localIndexSet0.index(i)[0];
for (int j=0; j<nDofs0; j++ ) {
- int col = basis0_.index(*it,j);
+ int col = localIndexSet0.index(j)[0];
hessian00[row][col] += localStiffness_->A00_[i][j];
}
for (int j=0; j<nDofs1; j++ ) {
- int col = basis1_.index(*it,j);
+ int col = localIndexSet1.index(j)[0];
hessian01[row][col] += localStiffness_->A01_[i][j];
}
}
for (int i=0; i<nDofs1; i++) {
- int row = basis1_.index(*it,i);
+ int row = localIndexSet1.index(i)[0];
for (int j=0; j<nDofs0; j++ ) {
- int col = basis0_.index(*it,j);
+ int col = localIndexSet0.index(j)[0];
hessian10[row][col] += localStiffness_->A10_[i][j];
}
for (int j=0; j<nDofs1; j++ ) {
- int col = basis1_.index(*it,j);
+ int col = localIndexSet1.index(j)[0];
hessian11[row][col] += localStiffness_->A11_[i][j];
}
}
// Add local gradient to global gradient
for (int i=0; i<nDofs0; i++)
- gradient0[basis0_.index(*it,i)] += localGradient0[i];
+ gradient0[localIndexSet0.index(i)[0]] += localGradient0[i];
for (int i=0; i<nDofs1; i++)
- gradient1[basis1_.index(*it,i)] += localGradient1[i];
+ gradient1[localIndexSet1.index(i)[0]] += localGradient1[i];
}
}
@@ -300,34 +315,46 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
{
double energy = 0;
- if (configuration0.size()!=basis0_.size())
- DUNE_THROW(Dune::Exception, "Configuration vector doesn't match the grid!");
+ if (configuration0.size()!=basis0_.indexSet().size())
+ DUNE_THROW(Dune::Exception, "Configuration vector 0 doesn't match the basis!");
+
+ if (configuration1.size()!=basis1_.indexSet().size())
+ DUNE_THROW(Dune::Exception, "Configuration vector 1 doesn't match the basis!");
- if (configuration1.size()!=basis1_.size())
- DUNE_THROW(Dune::Exception, "Configuration vector doesn't match the grid!");
+ // A view on the FE basis on a single element
+ typename Basis0::LocalView localView0(&basis0_);
+ typename Basis1::LocalView localView1(&basis1_);
+ auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+ auto localIndexSet1 = basis1_.indexSet().localIndexSet();
- ElementIterator it = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endIt = basis0_.getGridView().template end<0,Dune::Interior_Partition>();
+ ElementIterator it = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+ ElementIterator endIt = basis0_.gridView().template end<0,Dune::Interior_Partition>();
// Loop over all elements
for (; it!=endIt; ++it) {
+ // Bind the local FE basis view to the current element
+ localView0.bind(*it);
+ localView1.bind(*it);
+ localIndexSet0.bind(localView0);
+ localIndexSet1.bind(localView1);
+
// Number of degrees of freedom on this element
- size_t nDofs0 = basis0_.getLocalFiniteElement(*it).localBasis().size();
- size_t nDofs1 = basis1_.getLocalFiniteElement(*it).localBasis().size();
+ size_t nDofs0 = localView0.size();
+ size_t nDofs1 = localView1.size();
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (size_t i=0; i<nDofs0; i++)
- localConfiguration0[i] = configuration0[basis0_.index(*it,i)];
+ localConfiguration0[i] = configuration0[localIndexSet0.index(i)[0]];
for (size_t i=0; i<nDofs1; i++)
- localConfiguration1[i] = configuration1[basis1_.index(*it,i)];
+ localConfiguration1[i] = configuration1[localIndexSet1.index(i)[0]];
energy += localStiffness_->energy(*it,
- basis0_.getLocalFiniteElement(*it), localConfiguration0,
- basis1_.getLocalFiniteElement(*it), localConfiguration1);
+ localView0.tree().finiteElement(), localConfiguration0,
+ localView1.tree().finiteElement(), localConfiguration1);
}
diff --git a/dune/gfe/mixedlocalgeodesicfestiffness.hh b/dune/gfe/mixedlocalgeodesicfestiffness.hh
index 82a903570aca4624f3e9f3253f3cdbd89c8e4fea..bc0a827cb0fd8da3fe10631f800e4be5595f5f25 100644
--- a/dune/gfe/mixedlocalgeodesicfestiffness.hh
+++ b/dune/gfe/mixedlocalgeodesicfestiffness.hh
@@ -7,12 +7,17 @@
#include <dune/istl/matrix.hh>
-template<class GridView,
- class DeformationLocalFiniteElement, class DeformationTargetSpace,
- class OrientationLocalFiniteElement, class OrientationTargetSpace>
+template<class DeformationBasis, class DeformationTargetSpace,
+ class OrientationBasis, class OrientationTargetSpace>
class MixedLocalGeodesicFEStiffness
{
+ static_assert(std::is_same<typename DeformationBasis::GridView, typename OrientationBasis::GridView>::value,
+ "DeformationBasis and OrientationBasis must be designed on the same GridView!");
+
// grid types
+ typedef typename DeformationBasis::LocalView::Tree::FiniteElement DeformationLocalFiniteElement;
+ typedef typename OrientationBasis::LocalView::Tree::FiniteElement OrientationLocalFiniteElement;
+ typedef typename DeformationBasis::GridView GridView;
typedef typename GridView::Grid::ctype DT;
typedef typename DeformationTargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
diff --git a/dune/gfe/mixedlocalgfeadolcstiffness.hh b/dune/gfe/mixedlocalgfeadolcstiffness.hh
index a40212ca6781ee9a883a50a03120b1c1f531eff5..14d92b6adf05bcab90174becfdb181970ce61c29 100644
--- a/dune/gfe/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/mixedlocalgfeadolcstiffness.hh
@@ -18,15 +18,19 @@
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
-template<class GridView,
- class LocalFiniteElement0, class TargetSpace0,
- class LocalFiniteElement1, class TargetSpace1>
+template<class Basis0, class TargetSpace0,
+ class Basis1, class TargetSpace1>
class MixedLocalGFEADOLCStiffness
- : public MixedLocalGeodesicFEStiffness<GridView,
- LocalFiniteElement0,TargetSpace0,
- LocalFiniteElement1,TargetSpace1>
+ : public MixedLocalGeodesicFEStiffness<Basis0,TargetSpace0,
+ Basis1,TargetSpace1>
{
+ static_assert(std::is_same<typename Basis0::GridView, typename Basis1::GridView>::value,
+ "Basis0 and Basis1 must be designed on the same GridView!");
+
// grid types
+ typedef typename Basis0::GridView GridView;
+ typedef typename Basis0::LocalView::Tree::FiniteElement LocalFiniteElement0;
+ typedef typename Basis1::LocalView::Tree::FiniteElement LocalFiniteElement1;
typedef typename GridView::Grid::ctype DT;
typedef typename TargetSpace0::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -48,9 +52,8 @@ public:
enum { embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension };
enum { embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension };
- MixedLocalGFEADOLCStiffness(const MixedLocalGeodesicFEStiffness<GridView,
- LocalFiniteElement0, ATargetSpace0,
- LocalFiniteElement1, ATargetSpace1>* energy)
+ MixedLocalGFEADOLCStiffness(const MixedLocalGeodesicFEStiffness<Basis0, ATargetSpace0,
+ Basis1, ATargetSpace1>* energy)
: localEnergy_(energy)
{}
@@ -82,14 +85,14 @@ public:
std::vector<typename TargetSpace0::TangentVector>& localGradient0,
std::vector<typename TargetSpace1::TangentVector>& localGradient1);
- const MixedLocalGeodesicFEStiffness<GridView, LocalFiniteElement0, ATargetSpace0, LocalFiniteElement1, ATargetSpace1>* localEnergy_;
+ const MixedLocalGeodesicFEStiffness<Basis0, ATargetSpace0, Basis1, ATargetSpace1>* localEnergy_;
};
-template <class GridView, class LocalFiniteElement0, class TargetSpace0, class LocalFiniteElement1, class TargetSpace1>
-typename MixedLocalGFEADOLCStiffness<GridView, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::RT
-MixedLocalGFEADOLCStiffness<GridView, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::
+template <class Basis0, class TargetSpace0, class Basis1, class TargetSpace1>
+typename MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::RT
+MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::
energy(const Entity& element,
const LocalFiniteElement0& localFiniteElement0,
const std::vector<TargetSpace0>& localConfiguration0,
@@ -198,8 +201,8 @@ assembleGradient(const Entity& element,
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
-template <class GridType, class LocalFiniteElement0, class TargetSpace0, class LocalFiniteElement1, class TargetSpace1>
-void MixedLocalGFEADOLCStiffness<GridType, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::
+template <class Basis0, class TargetSpace0, class Basis1, class TargetSpace1>
+void MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement0& localFiniteElement0,
const std::vector<TargetSpace0>& localConfiguration0,
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 32cf12dfd46c878c5a62b697cb52a22bde327b25..b334c8a4080b161aa4d723b6939b7c8e76bdda56 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -7,17 +7,16 @@
#include <dune/istl/io.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
// Using a monotone multigrid as the inner solver
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
-#if defined THIRD_ORDER || defined SECOND_ORDER
-#include <dune/gfe/pktop1mgtransfer.hh>
-#endif
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include "maxnormtrustregion.hh"
@@ -135,11 +134,16 @@ setup(const GridType& grid,
mmgStep1->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType1>();
mmgStep1->verbosity_ = Solver::FULL;
-#if 0
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
+ typedef DuneFunctionsBasis<Basis0> FufemBasis0;
+ FufemBasis0 basis0(grid.leafGridView());
+
+ typedef DuneFunctionsBasis<Basis1> FufemBasis1;
+ FufemBasis1 basis1(grid.leafGridView());
+#if 0
BasisType basis(grid.leafGridView());
OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
@@ -189,47 +193,27 @@ setup(const GridType& grid,
mmgStep0->mgTransfer_.resize(numLevels-1);
mmgStep1->mgTransfer_.resize(numLevels-1);
- if (assembler->basis0_.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
+ if (basis0.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
{
if (numLevels>1) {
+ typedef typename TruncatedCompressedMGTransfer<CorrectionType0>::TransferOperatorType TransferOperatorType;
P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
- PKtoP1MGTransfer<CorrectionType0>* topTransferOp0 = new PKtoP1MGTransfer<CorrectionType0>;
- topTransferOp0->setup(assembler->basis0_,p1Basis);
- #if 1
- mmgStep0->mgTransfer_.back() = topTransferOp0;
- #else
- // If we are on more than 1 processors, join all local transfer matrices on rank 0,
- // and construct a single global transfer operator there.
- typedef GlobalUniqueIndex<typename GridType::LeafGridView, gridDim> LeafP1GUIndex;
- LeafP1GUIndex p1Index(grid_->leafGridView());
-
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- MatrixCommunicator<GUIndex, TransferOperatorType, LeafP1GUIndex> matrixComm(*guIndex_, p1Index, 0);
-
- mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>
- (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix())));
- #endif
- for (int i=0; i<mmgStep0->mgTransfer_.size()-1; i++){
+ TransferOperatorType pkToP1TransferMatrix;
+ assembleBasisInterpolationMatrix<TransferOperatorType,
+ P1NodalBasis<typename GridType::LeafGridView,double>,
+ FufemBasis0>(pkToP1TransferMatrix,p1Basis,assembler->basis0_);
+
+ mmgStep0->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType0>;
+ Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+ dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType0>*>(mmgStep0->mgTransfer_.back())->setMatrix(topTransferOperator);
+
+ for (size_t i=0; i<mmgStep0->mgTransfer_.size()-1; i++){
// Construct the local multigrid transfer matrix
TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
newTransferOp0->setup(*grid_,i+1,i+2);
- #if 1
mmgStep0->mgTransfer_[i] = newTransferOp0;
- #else
- // If we are on more than 1 processors, join all local transfer matrices on rank 0,
- // and construct a single global transfer operator there.
- typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
- LevelGUIndex fineGUIndex(grid_->levelGridView(i+2));
- LevelGUIndex coarseGUIndex(grid_->levelGridView(i+1));
-
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
-
- mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
- (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
- #endif
}
}
@@ -244,35 +228,27 @@ setup(const GridType& grid,
TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
newTransferOp0->setup(*grid_,i,i+1);
- #if 1
mmgStep0->mgTransfer_[i] = newTransferOp0;
- #else
- // If we are on more than 1 processors, join all local transfer matrices on rank 0,
- // and construct a single global transfer operator there.
- typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
- LevelGUIndex fineGUIndex(grid_->levelGridView(i+1));
- LevelGUIndex coarseGUIndex(grid_->levelGridView(i));
-
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
-
- mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
- (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
- #endif
}
}
- if (assembler->basis1_.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
+ if (basis1.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
{
if (numLevels>1) {
+ typedef typename TruncatedCompressedMGTransfer<CorrectionType1>::TransferOperatorType TransferOperatorType;
P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
- PKtoP1MGTransfer<CorrectionType1>* topTransferOp1 = new PKtoP1MGTransfer<CorrectionType1>;
- topTransferOp1->setup(assembler->basis1_,p1Basis);
- mmgStep1->mgTransfer_.back() = topTransferOp1;
+ TransferOperatorType pkToP1TransferMatrix;
+ assembleBasisInterpolationMatrix<TransferOperatorType,
+ P1NodalBasis<typename GridType::LeafGridView,double>,
+ FufemBasis1>(pkToP1TransferMatrix,p1Basis,assembler->basis1_);
+
+ mmgStep0->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType1>;
+ Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+ dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType1>*>(mmgStep1->mgTransfer_.back())->setMatrix(topTransferOperator);
- for (int i=0; i<mmgStep1->mgTransfer_.size()-1; i++){
+ for (size_t i=0; i<mmgStep1->mgTransfer_.size()-1; i++){
// Construct the local multigrid transfer matrix
TruncatedCompressedMGTransfer<CorrectionType1>* newTransferOp1 = new TruncatedCompressedMGTransfer<CorrectionType1>;
newTransferOp1->setup(*grid_,i+1,i+2);
@@ -303,8 +279,8 @@ setup(const GridType& grid,
#if 0
hasObstacle0_.resize(guIndex_->nGlobalEntity(), true);
#else
- hasObstacle0_.resize(assembler->basis0_.size(), true);
- hasObstacle1_.resize(assembler->basis1_.size(), true);
+ hasObstacle0_.resize(assembler->basis0_.indexSet().size(), true);
+ hasObstacle1_.resize(assembler->basis1_.indexSet().size(), true);
#endif
mmgStep0->hasObstacle_ = &hasObstacle0_;
mmgStep1->hasObstacle_ = &hasObstacle1_;
@@ -323,15 +299,9 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc,argv);
int rank = grid_->comm().rank();
- MonotoneMGStep<MatrixType00,CorrectionType0>* mgStep0 = nullptr;
-
- // if the inner solver is a monotone multigrid set up a max-norm trust-region
- if (dynamic_cast<LoopSolver<CorrectionType0>*>(innerSolver_.get()))
- mgStep0 = dynamic_cast<MonotoneMGStep<MatrixType00,CorrectionType0>*>(dynamic_cast<LoopSolver<CorrectionType0>*>(innerSolver_.get())->iterationStep_);
-
// \todo Use global index set instead of basis for parallel computations
- MaxNormTrustRegion<blocksize0> trustRegion0(assembler_->basis0_.size(), initialTrustRegionRadius_);
- MaxNormTrustRegion<blocksize1> trustRegion1(assembler_->basis1_.size(), initialTrustRegionRadius_);
+ MaxNormTrustRegion<blocksize0> trustRegion0(assembler_->basis0_.indexSet().size(), initialTrustRegionRadius_);
+ MaxNormTrustRegion<blocksize1> trustRegion1(assembler_->basis1_.indexSet().size(), initialTrustRegionRadius_);
std::vector<BoxConstraint<field_type,blocksize0> > trustRegionObstacles0;
std::vector<BoxConstraint<field_type,blocksize1> > trustRegionObstacles1;
@@ -607,10 +577,12 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
}
// Output each iterate, to better understand what the algorithm does
+ DuneFunctionsBasis<Basis0> fufemBasis0(assembler_->basis0_);
+ DuneFunctionsBasis<Basis1> fufemBasis1(assembler_->basis1_);
std::stringstream iAsAscii;
iAsAscii << i+1;
- CosseratVTKWriter<GridType>::template writeMixed<Basis0,Basis1>(assembler_->basis0_,x0_,
- assembler_->basis1_,x1_,
+ CosseratVTKWriter<GridType>::template writeMixed<DuneFunctionsBasis<Basis0>, DuneFunctionsBasis<Basis1> >(fufemBasis0,x0_,
+ fufemBasis1,x1_,
"mixed-cosserat_iterate_" + iAsAscii.str());
if (rank==0)
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index d3c64fc2fb034aa2a012b6494b2e361c9b79a917..df3186ffb8bd759ce10b105888be5304f520c793 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -11,25 +11,32 @@
-template <class GridView>
-void RodAssembler<GridView,3>::
+template <class Basis>
+void RodAssembler<Basis,3>::
assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
using namespace Dune;
- if (sol.size()!=this->basis_.size())
+ if (sol.size()!=this->basis_.indexSet().size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
- ElementIterator it = this->basis_.getGridView().template begin<0>();
- ElementIterator endIt = this->basis_.getGridView().template end<0>();
+ // A view on the FE basis on a single element
+ typename Basis::LocalView localView(&this->basis_);
+ auto localIndexSet = this->basis_.indexSet().localIndexSet();
+
+ ElementIterator it = this->basis_.gridView().template begin<0>();
+ ElementIterator endIt = this->basis_.gridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
+ localView.bind(*it);
+ localIndexSet.bind(localView);
+
// A 1d grid has two vertices
static const int nDofs = 2;
@@ -37,19 +44,19 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
- localSolution[i] = sol[this->basis_.index(*it,i)];
+ localSolution[i] = sol[localIndexSet.index(i)[0]];
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
this->localStiffness_->assembleGradient(*it,
- this->basis_.getLocalFiniteElement(*it),
+ localView.tree().finiteElement(),
localSolution,
localGradient);
// Add to global gradient
for (int i=0; i<nDofs; i++)
- grad[this->basis_.index(*it,i)] += localGradient[i];
+ grad[localIndexSet.index(i)[0]] += localGradient[i];
}
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index a40dcb99e89c6280e732b0a7c04ca5b2fca68e8a..40fadb719b3b62b5a1ce215fbc6dbe412364cb22 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -14,7 +14,7 @@
/** \brief The FEM operator for an extensible, shearable rod in 3d
*/
-template <class GridView, int spaceDim>
+template <class Basis, int spaceDim>
class RodAssembler
{
static_assert(spaceDim==2 || spaceDim==3,
@@ -23,12 +23,11 @@ class RodAssembler
/** \brief The FEM operator for an extensible, shearable rod in 3d
*/
-template <class GridView>
-class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >
+template <class Basis>
+class RodAssembler<Basis,3> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >
{
+ typedef typename Basis::GridView GridView;
- //typedef typename GridType::template Codim<0>::Entity EntityType;
- //typedef typename GridType::template Codim<0>::EntityPointer EntityPointer;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
//! Dimension of the grid. This needs to be one!
@@ -44,18 +43,18 @@ class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridVie
public:
//! ???
- RodAssembler(const GridView &gridView,
+ RodAssembler(const Basis& basis,
RodLocalStiffness<GridView,double>* localStiffness)
- : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >(gridView,localStiffness)
+ : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
{
- std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(gridDim));
+ std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.indexSet().size());
- typename GridView::template Codim<gridDim>::Iterator it = gridView.template begin<gridDim>();
- typename GridView::template Codim<gridDim>::Iterator endIt = gridView.template end<gridDim>();
+ auto it = basis.gridView().template begin<gridDim>();
+ auto endIt = basis.gridView().template end<gridDim>();
for (; it != endIt; ++it) {
- int idx = gridView.indexSet().index(*it);
+ int idx = basis.gridView().indexSet().index(*it);
referenceConfiguration[idx].r[0] = 0;
referenceConfiguration[idx].r[1] = 0;
@@ -91,10 +90,11 @@ public:
/** \brief The FEM operator for a 2D extensible, shearable rod
*/
-template <class GridView>
-class RodAssembler<GridView,2> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >
+template <class Basis>
+class RodAssembler<Basis,2> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >
{
+ typedef typename Basis::GridView GridView;
typedef typename GridView::template Codim<0>::Entity EntityType;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -118,7 +118,7 @@ public:
//! ???
RodAssembler(const GridView &gridView)
- : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >(gridView,NULL)
+ : GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >(gridView,NULL)
{
B = 1;
A1 = 1;
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index 3b03fa8057a90e073177025f1d17eec2bbb737c6..e191817aabd602df8418f56ec96dce3308ae5059 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -4,16 +4,15 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
#include <dune/geometry/quadraturerules.hh>
-#include <dune/localfunctions/lagrange/p1.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include "localgeodesicfestiffness.hh"
#include "rigidbodymotion.hh"
template<class GridView, class RT>
class RodLocalStiffness
- : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<RT,3> >
+ : public LocalGeodesicFEStiffness<Dune::Functions::PQKNodalBasis<GridView,1>, RigidBodyMotion<RT,3> >
{
typedef RigidBodyMotion<RT,3> TargetSpace;
@@ -102,7 +101,7 @@ public:
const Dune::array<RigidBodyMotion<RT,3>, dim+1>& localSolution) const;
virtual RT energy (const Entity& e,
- const typename P1NodalBasis<GridView>::LocalFiniteElement& localFiniteElement,
+ const typename Dune::Functions::PQKNodalBasis<GridView,1>::LocalView::Tree::FiniteElement& localFiniteElement,
const std::vector<RigidBodyMotion<RT,3> >& localSolution) const
{
assert(localSolution.size()==2);