diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 1a5005e66a4f6b7a95bdf9014ea36b4cedcb479f..e8881dcacd920ebd2b8d48627d0b914ca29655a9 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -662,7 +662,7 @@ assembleMatrixFD(const std::vector<Configuration>& sol,
Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<GridType,double> localStiffness(K, A);
+ RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
// //////////////////////////////////////////////////////////////////
// Store pointers to all elements so we can access them by index
@@ -904,7 +904,7 @@ assembleGradient(const std::vector<Configuration>& sol,
Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<GridType,double> localStiffness(K, A);
+ RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
grad.resize(sol.size());
grad = 0;
@@ -963,7 +963,7 @@ computeEnergy(const std::vector<Configuration>& sol) const
Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<GridType,double> localStiffness(K, A);
+ RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
Dune::array<Configuration,2> localReferenceConfiguration;
Dune::array<Configuration,2> localSolution;
@@ -1008,7 +1008,7 @@ getStrain(const std::vector<Configuration>& sol,
Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<GridType,double> localStiffness(K, A);
+ RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
// Strain defined on each element
strain.resize(indexSet.size(0));
@@ -1128,7 +1128,7 @@ getResultantForce(const BoundaryPatch<GridType>& boundary,
Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<GridType,double> localStiffness(K, A);
+ RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
double pos = nIt->intersectionSelfLocal()[0];
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 80beadf3a228c944800899149df36ae51967e583..3fd4cbc5cea62bf5786969d0687d3b597a397284 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -10,17 +10,17 @@
#include <dune/ag-common/boundarypatch.hh>
#include "configuration.hh"
-template<class GridType, class RT>
+template<class GridView, class RT>
class RodLocalStiffness
- : public Dune::LocalStiffness<RodLocalStiffness<GridType,RT>,GridType,RT,6>
+ : public Dune::LocalStiffness<GridView,RT,6>
{
// grid types
- typedef typename GridType::ctype DT;
- typedef typename GridType::template Codim<0>::Entity Entity;
- typedef typename GridType::template Codim<0>::EntityPointer EntityPointer;
+ typedef typename GridView::Grid::ctype DT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+ typedef typename GridView::template Codim<0>::EntityPointer EntityPointer;
// some other sizes
- enum {dim=GridType::dimension};
+ enum {dim=GridView::dimension};
// Quadrature order used for the extension and shear energy
enum {shearQuadOrder = 2};
@@ -37,7 +37,7 @@ public:
// to allocate the correct size of (dense) blocks with a FieldMatrix
enum {m=blocksize};
- enum {SIZE = Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE};
+ enum {SIZE = Dune::LocalStiffness<GridView,RT,m>::SIZE};
// types for matrics, vectors and boundary conditions
typedef Dune::FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
@@ -59,7 +59,7 @@ public:
// For the time being: all boundary conditions are homogeneous Neumann
// This means no boundary condition handling is done at all
- for (int i=0; i<Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE; i++)
+ for (int i=0; i<SIZE; i++)
for (size_t j=0; j<this->bctype[i].size(); j++)
this->bctype[i][j] = Dune::BoundaryConditions::neumann;
}
@@ -75,7 +75,7 @@ public:
// For the time being: all boundary conditions are homogeneous Neumann
// This means no boundary condition handling is done at all
- for (int i=0; i<Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE; i++)
+ for (int i=0; i<SIZE; i++)
for (size_t j=0; j<this->bctype[i].size(); j++)
this->bctype[i][j] = Dune::BoundaryConditions::neumann;
}
@@ -91,10 +91,23 @@ public:
\param[in] localSolution Current local solution, because this is a nonlinear assembler
@param[in] k order of Lagrange basis
*/
- template<typename TypeTag>
void assemble (const Entity& e,
- const Dune::BlockVector<Dune::FieldVector<double, dim> >& localSolution,
- int k=1);
+ const Dune::BlockVector<Dune::FieldVector<double, 6> >& localSolution,
+ int k=1)
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
+
+ /** \todo Remove this once this methods is not in base class LocalStiffness anymore */
+ void assemble (const Entity& e, int k=1)
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
+
+ void assembleBoundaryCondition (const Entity& e, int k=1)
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
RT energy (const Entity& e,
diff --git a/src/rodsolver.cc b/src/rodsolver.cc
index 5ecd5f8420640ef206d95043a9c3b2944f81df8f..88d7d1c8fadad53ad3457614d56097e24a42e9e7 100644
--- a/src/rodsolver.cc
+++ b/src/rodsolver.cc
@@ -49,6 +49,7 @@ template <class GridType>
void RodSolver<GridType>::setup(const GridType& grid,
const RodAssembler<GridType>* rodAssembler,
const SolutionType& x,
+ const BlockBitField<blocksize>& dirichletNodes,
double tolerance,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
@@ -80,26 +81,11 @@ void RodSolver<GridType>::setup(const GridType& grid,
int numLevels = grid_->maxLevel()+1;
- // ////////////////////////////////////////////
- // Construct array with the Dirichlet nodes
- // ////////////////////////////////////////////
- dirichletNodes_.resize(numLevels);
- for (int i=0; i<numLevels; i++) {
-
- dirichletNodes_[i].resize( blocksize * grid.size(i,1), false );
-
- for (int j=0; j<blocksize; j++) {
- dirichletNodes_[i][j] = true;
- dirichletNodes_[i][dirichletNodes_[i].size()-1-j] = true;
- }
-
- }
-
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
- // First create a gauss-seidel base solver
+ // First create a Gauss-seidel base solver
TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
@@ -118,7 +104,7 @@ void RodSolver<GridType>::setup(const GridType& grid,
MonotoneMGStep<MatrixType, CorrectionType>* mmgStep = new MonotoneMGStep<MatrixType, CorrectionType>(numLevels);
mmgStep->setMGType(mu_, nu1_, nu2_);
- mmgStep->dirichletNodes_ = &dirichletNodes_;
+ mmgStep->ignoreNodes = &dirichletNodes;
mmgStep->basesolver_ = baseSolver;
mmgStep->presmoother_ = presmoother;
mmgStep->postsmoother_ = postsmoother;
@@ -131,7 +117,7 @@ void RodSolver<GridType>::setup(const GridType& grid,
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
LeafP1Function<GridType,double> u(grid),f(grid);
- LaplaceLocalStiffness<GridType,double> laplaceStiffness;
+ LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
LeafP1OperatorAssembler<GridType,double,1>* A = new LeafP1OperatorAssembler<GridType,double,1>(grid);
A->assemble(laplaceStiffness,u,f);
diff --git a/src/rodsolver.hh b/src/rodsolver.hh
index 4850ac68f838df0a24ed65079df1688c97f3eeba..5fec9522ffbbf8bc9d0074f5e830166096ebfd9f 100644
--- a/src/rodsolver.hh
+++ b/src/rodsolver.hh
@@ -9,6 +9,7 @@
#include <dune/ag-common/boxconstraint.hh>
#include <dune/ag-common/norms/h1seminorm.hh>
#include <dune/ag-common/iterativesolver.hh>
+#include <dune/ag-common/blockbitfield.hh>
#include "rodassembler.hh"
@@ -40,6 +41,7 @@ public:
void setup(const GridType& grid,
const RodAssembler<GridType>* rodAssembler,
const SolutionType& x,
+ const BlockBitField<blocksize>& dirichletNodes,
double tolerance,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
@@ -110,9 +112,6 @@ protected:
expects them :-( */
std::vector<Dune::BitField> hasObstacle_;
- /** \brief The Dirichlet nodes on all levels */
- std::vector<Dune::BitField> dirichletNodes_;
-
/** \brief The norm used to measure multigrid convergence */
H1SemiNorm<CorrectionType>* h1SemiNorm_;