diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index 8e852b731c76c86d38297c2a6302d530ef427b3f..c6d488552ee35c6067784775cdeb3f949e4bc959 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -3,13 +3,10 @@
#include <dune/grid/onedgrid.hh>
#include <dune/grid/uggrid.hh>
-#include <dune/disc/elasticity/linearelasticityassembler.hh>
-#include <dune/disc/operators/p1operator.hh>
#include <dune/istl/io.hh>
#include <dune/grid/io/file/amirameshreader.hh>
#include <dune/grid/io/file/amirameshwriter.hh>
-
#include <dune/common/bitsetvector.hh>
#include <dune/common/configparser.hh>
@@ -27,6 +24,10 @@
#include <dune/ag-common/neumannassembler.hh>
#include <dune/ag-common/computestress.hh>
+#include <dune/ag-common/functionspacebases/q1nodalbasis.hh>
+#include <dune/ag-common/assemblers/operatorassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+
#include "src/quaternion.hh"
#include "src/rodassembler.hh"
#include "src/rigidbodymotion.hh"
@@ -238,10 +239,15 @@ int main (int argc, char *argv[]) try
// //////////////////////////////////////////
// Assemble 3d linear elasticity problem
// //////////////////////////////////////////
- LeafP1Function<GridType,double,dim> u(grid),f(grid);
- LinearElasticityLocalStiffness<GridType::LeafGridView,double> lstiff(E, nu);
- LeafP1OperatorAssembler<GridType,double,dim> hessian3d(grid);
- hessian3d.assemble(lstiff,u,f);
+
+ typedef Q1NodalBasis<GridType::LeafGridView,double> FEBasis;
+ FEBasis basis(grid.leafView());
+ OperatorAssembler<FEBasis,FEBasis> assembler(basis, basis);
+
+ StVenantKirchhoffAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> localAssembler(E, nu);
+ MatrixType stiffnessMatrix3d;
+
+ assembler.assemble(localAssembler, stiffnessMatrix3d);
// ////////////////////////////////////////////////////////////
// Create solution and rhs vectors
@@ -318,7 +324,7 @@ int main (int argc, char *argv[]) try
// Make pre and postsmoothers
BlockGSStep<MatrixType, VectorType> presmoother, postsmoother;
- MultigridStep<MatrixType, VectorType> multigridStep(*hessian3d, x3d, rhs3d, 1);
+ MultigridStep<MatrixType, VectorType> multigridStep(stiffnessMatrix3d, x3d, rhs3d, 1);
multigridStep.setMGType(mu, nu1, nu2);
multigridStep.ignoreNodes_ = &dirichletNodes.back();
@@ -419,7 +425,7 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////////////////////////
// Solve the Neumann problem for the 3d body
// ///////////////////////////////////////////////////////////
- multigridStep.setProblem(*hessian3d, x3d, rhs3d, grid.maxLevel()+1);
+ multigridStep.setProblem(stiffnessMatrix3d, x3d, rhs3d, grid.maxLevel()+1);
solver.preprocess();
multigridStep.preprocess();
@@ -500,9 +506,9 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////
// the 3d body
- double oldNorm = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, *hessian3d);
+ double oldNorm = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
oldSolution3d -= x3d;
- double normOfCorrection = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, *hessian3d);
+ double normOfCorrection = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
double max3dRelCorrection = 0;
for (size_t j=0; j<x3d.size(); j++)
@@ -581,7 +587,7 @@ int main (int argc, char *argv[]) try
// This should really be exactSol-initialSol, but we're starting
// from zero anyways
- oldError += EnergyNorm<MatrixType,VectorType>::normSquared(exactSol3d, *hessian3d);
+ oldError += EnergyNorm<MatrixType,VectorType>::normSquared(exactSol3d, stiffnessMatrix3d);
// Error of the initial rod iterate
RodDifferenceType rodDifference = computeGeodesicDifference(initialIterateRod, exactSolRod);
@@ -646,7 +652,7 @@ int main (int argc, char *argv[]) try
RodDifferenceType rodDifference = computeGeodesicDifference(exactSolRod, intermediateSolRod);
- error = std::sqrt(EnergyNorm<MatrixType,VectorType>::normSquared(solBackup0, *hessian3d)
+ error = std::sqrt(EnergyNorm<MatrixType,VectorType>::normSquared(solBackup0, stiffnessMatrix3d)
+
EnergyNorm<BCRSMatrix<FieldMatrix<double,6,6> >,RodDifferenceType>::normSquared(rodDifference, hessianRod));
@@ -712,7 +718,7 @@ int main (int argc, char *argv[]) try
amiraMeshWriter.addVertexData(x3d, grid.leafView());
BlockVector<FieldVector<double,1> > stress;
- Stress<GridType,dim>::getStress(grid, x3d, stress, E, nu);
+ Stress<GridType>::getStress(grid, x3d, stress, E, nu);
amiraMeshWriter.addVertexData(stress, grid.leafView());
amiraMeshWriter.write(resultPath + "grid.result");
diff --git a/rod-eoc.cc b/rod-eoc.cc
index 67b5f2eef8be772b0d6e9a59dd455845c0419c28..1baee43f179811130ade47b2ecfa0573d63b24f7 100644
--- a/rod-eoc.cc
+++ b/rod-eoc.cc
@@ -7,8 +7,10 @@
#include <dune/istl/io.hh>
-#include <dune/disc/miscoperators/massmatrix.hh>
-#include <dune/disc/miscoperators/laplace.hh>
+#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
+#include <dune/ag-common/assemblers/operatorassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/massassembler.hh>
#include <dune-solvers/solvers/iterativesolver.hh>
#include <dune-solvers/norms/energynorm.hh>
@@ -179,15 +181,18 @@ int main (int argc, char *argv[]) try
// Compute mass matrix and laplace matrix to emulate L2 and H1 norms
// //////////////////////////////////////////////////////////////////////
- Dune::LeafP1Function<GridType,double> u(referenceGrid),f(referenceGrid);
+ typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
+ FEBasis basis(referenceGrid.leafView());
+ OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
- Dune::MassMatrixLocalStiffness<GridType::LeafGridView,double,1> massMatrixStiffness;
- Dune::LeafP1OperatorAssembler<GridType,double,1> massMatrix(referenceGrid);
- massMatrix.assemble(massMatrixStiffness,u,f);
+ LaplaceAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> laplaceLocalAssembler;
+ MassAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> massMatrixLocalAssembler;
- Dune::LaplaceLocalStiffness<GridType::LeafGridView,double> laplaceStiffness;
- Dune::LeafP1OperatorAssembler<GridType,double,1> laplace(referenceGrid);
- laplace.assemble(laplaceStiffness,u,f);
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+ ScalarMatrixType laplace, massMatrix;
+
+ operatorAssembler.assemble(laplaceLocalAssembler, laplace);
+ operatorAssembler.assemble(massMatrixLocalAssembler, massMatrix);
// ///////////////////////////////////////////////////////////
// Compute on all coarser levels, and compare
@@ -214,8 +219,8 @@ int main (int argc, char *argv[]) try
BlockVector<TargetSpace::TangentVector> difference = computeGeodesicDifference(solution,referenceSolution);
- H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(*laplace);
- H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(*massMatrix);
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(laplace);
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(massMatrix);
// Compute max-norm difference
std::cout << "Level: " << i-1
diff --git a/src/riemanniantrsolver.cc b/src/riemanniantrsolver.cc
index 4647c6cb79d0f948ca518dde75b28745e940a0de..bbc9259590089545eb071a3c5aee6c591c378fc4 100644
--- a/src/riemanniantrsolver.cc
+++ b/src/riemanniantrsolver.cc
@@ -2,10 +2,10 @@
#include <dune/istl/io.hh>
-#include <dune/disc/functions/p1function.hh>
-#include <dune/disc/operators/p1operator.hh>
-#include <dune/disc/miscoperators/laplace.hh>
-#include <dune/disc/miscoperators/massmatrix.hh>
+#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
+#include <dune/ag-common/assemblers/operatorassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/massassembler.hh>
// For using a monotone multigrid as the inner solver
#include <dune-solvers/iterationsteps/trustregiongsstep.hh>
@@ -89,17 +89,20 @@ setup(const GridType& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
- Dune::LeafP1Function<GridType,double> u(grid),f(grid);
- Dune::LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
- Dune::LeafP1OperatorAssembler<GridType,double,1>* A = new Dune::LeafP1OperatorAssembler<GridType,double,1>(grid);
- A->assemble(laplaceStiffness,u,f);
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
+ FEBasis basis(grid.leafView());
+ OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
- typedef typename Dune::LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
+ LaplaceAssembler<GridType, typename FEBasis::LocalFiniteElement, typename FEBasis::LocalFiniteElement> laplaceStiffness;
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+ ScalarMatrixType* A = new ScalarMatrixType;
+
+ operatorAssembler.assemble(laplaceStiffness, *A);
if (h1SemiNorm_)
delete h1SemiNorm_;
- h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A);
+ h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
innerSolver_ = new ::LoopSolver<CorrectionType>(mmgStep,
innerIterations_,
@@ -183,25 +186,30 @@ setupTCG(const GridType& grid,
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// This is used to measure convergence of the inner solver
// //////////////////////////////////////////////////////////////////////////////////////
- Dune::LeafP1Function<GridType,double> u(grid),f(grid);
- Dune::LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
- Dune::LeafP1OperatorAssembler<GridType,double,1>* A = new Dune::LeafP1OperatorAssembler<GridType,double,1>(grid);
- A->assemble(laplaceStiffness,u,f);
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
+ FEBasis basis(grid.leafView());
+ OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
+
+ LaplaceAssembler<GridType, typename FEBasis::LocalFiniteElement, typename FEBasis::LocalFiniteElement> laplaceStiffness;
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+ ScalarMatrixType* A = new ScalarMatrixType;
+
+ operatorAssembler.assemble(laplaceStiffness, *A);
if (h1SemiNorm_)
delete h1SemiNorm_;
- h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A);
+ h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a mass matrix to create a norm that's equivalent to the L2-norm
// This is used to to define the trust region
// //////////////////////////////////////////////////////////////////////////////////////
- Dune::LeafP1Function<GridType,double,blocksize> uvv(grid),fvv(grid);
- Dune::MassMatrixLocalStiffness<typename GridType::LeafGridView,double,blocksize> massMatrixStiffness;
- Dune::LeafP1OperatorAssembler<GridType,double,blocksize>* B = new Dune::LeafP1OperatorAssembler<GridType,double,blocksize>(grid);
- B->assemble(massMatrixStiffness,uvv,fvv);
+ MassAssembler<GridType, typename FEBasis::LocalFiniteElement, typename FEBasis::LocalFiniteElement> massMatrixStiffness;
+ MatrixType* B = new MatrixType;
+
+ operatorAssembler.assemble(massMatrixStiffness, *B);
// /////////////////////////////////////////////////////////
// Scale rotational components
@@ -209,11 +217,11 @@ setupTCG(const GridType& grid,
#warning RigidBody-specific stuff hardwired into the RiemannianTRSolver!
int alpha = 1;
- for (int i=0; i<(**B).N(); i++) {
+ for (int i=0; i<B->N(); i++) {
// make matrix row an identity row
- typename MatrixType::row_type::Iterator cIt = (**B)[i].begin();
- typename MatrixType::row_type::Iterator cEndIt = (**B)[i].end();
+ typename MatrixType::row_type::Iterator cIt = (*B)[i].begin();
+ typename MatrixType::row_type::Iterator cEndIt = (*B)[i].end();
for (; cIt!=cEndIt; ++cIt)
for (int j=0; j<3; j++)
@@ -229,7 +237,7 @@ setupTCG(const GridType& grid,
innerIterations_,
innerTolerance_,
h1SemiNorm_,
- &**B, // the norm of the trust region
+ B, // the norm of the trust region
Solver::FULL);
// Write all intermediate solutions, if requested