a very crude test for the MonotoneMGStep class

[[Imported from SVN: r7396]]

dune/solvers/test/Makefile.am

@@ -2,6 +2,7 @@
 # list of tests to run
 TESTS = genericvectortoolstest \
         lowrankoperatortest \
+        mmgtest \
         multigridtest \
         nulloperatortest \
         obstacletnnmgtest \
@@ -34,6 +35,11 @@ lowrankoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
 lowrankoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
 lowrankoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
 
+mmgtest_SOURCES = mmgtest.cc
+mmgtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) $(IPOPT_CPPFLAGS)
+mmgtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) $(IPOPT_LIBS)
+mmgtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) $(IPOPT_LIBS)
+
 multigridtest_SOURCES = multigridtest.cc
 multigridtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
 multigridtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)

dune/solvers/test/mmgtest.cc

+#include <config.h>
+
+#include <iostream>
+#include <sstream>
+
+// dune-common includes
+#include <dune/common/bitsetvector.hh>
+
+// dune-istl includes
+#include <dune/istl/bcrsmatrix.hh>
+
+// dune-grid includes
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+// dune-solver includes
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
+
+
+#include "common.hh"
+
+template <class DomainType, class RangeType>
+class ConstantFunction :
+    public Dune::VirtualFunction<DomainType, RangeType>
+{
+    public:
+        ConstantFunction(double c):
+            c_(c)
+        {}
+
+        void evaluate(const DomainType& x, RangeType& y) const
+        {
+            y = c_;
+        }
+
+    private:
+        double c_;
+};
+
+
+template<class GridType, class MatrixType, class VectorType, class BitVector>
+void solveObstacleProblemByMMGSolver(const GridType& grid, const MatrixType& mat, VectorType& x, const VectorType& rhs, const BitVector& ignore, int maxIterations=100, double tolerance=1.0e-10)
+{
+    typedef VectorType Vector;
+    typedef MatrixType Matrix;
+    typedef EnergyNorm<Matrix, Vector> Norm;
+    typedef QuadraticIPOptSolver<MatrixType,VectorType> Solver;
+
+    // create double obstacle constraints
+#if 0
+    const int blockSize = VectorType::block_type::dimension;
+
+    BoxConstraintVector boxConstraints(rhs.size());
+    for (int i = 0; i < boxConstraints.size(); ++i)
+    {
+        for (int j = 0; j < blockSize; ++j)
+        {
+            boxConstraints[i].lower(j) = -1;
+            boxConstraints[i].upper(j) = +1;
+        }
+    }
+#endif
+
+    /////////////////////////////////////////////////////////////
+    //  Create a monotone multigrid solver
+    /////////////////////////////////////////////////////////////
+    
+    // First create a Gauss-seidel base solver
+    ProjectedBlockGSStep<MatrixType, VectorType>* baseSolverStep = new ProjectedBlockGSStep<MatrixType, VectorType>;
+
+    EnergyNorm<MatrixType, VectorType>* baseEnergyNorm = new EnergyNorm<MatrixType, VectorType>(*baseSolverStep);
+
+    ::LoopSolver<VectorType>* baseSolver = new ::LoopSolver<VectorType>(baseSolverStep,
+                                                                        100,   // base iterations
+                                                                        1e-8,  // tolerance
+                                                                        baseEnergyNorm,
+                                                                        Solver::QUIET);
+    
+    // Make pre and postsmoothers
+    ProjectedBlockGSStep<MatrixType, VectorType>* presmoother  = new ProjectedBlockGSStep<MatrixType, VectorType>;
+    ProjectedBlockGSStep<MatrixType, VectorType>* postsmoother = new ProjectedBlockGSStep<MatrixType, VectorType>;
+
+    MonotoneMGStep<MatrixType, VectorType>* mmgStep = new MonotoneMGStep<MatrixType, VectorType>();
+
+    mmgStep->setMGType(3, 3, 1);
+    mmgStep->ignoreNodes_       = &ignore;
+    mmgStep->basesolver_        = baseSolver;
+    mmgStep->setSmoother(presmoother, postsmoother);
+    mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<VectorType>();
+    mmgStep->hasObstacle_       = NULL;
+    mmgStep->verbosity_         = Solver::QUIET;
+
+    EnergyNorm<MatrixType,VectorType>* h1SemiNorm = new EnergyNorm<MatrixType,VectorType>(mat);
+
+    LoopSolver<VectorType> solver (mmgStep,
+                                   100,    // iterations
+                                   1e-6,   // tolerance
+                                   h1SemiNorm,
+                                   Solver::REDUCED);
+    // solve problem
+    solver.solve();
+}
+
+
+
+
+
+template <class GridType>
+bool checkWithGrid(const GridType& grid, const std::string fileName="")
+{
+    bool passed = true;
+
+    static const int dim = GridType::dimension;
+
+    typedef typename Dune::FieldMatrix<double, 1, 1> MatrixBlock;
+    typedef typename Dune::FieldVector<double, 1> VectorBlock;
+    typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
+    typedef typename Dune::BlockVector<VectorBlock> Vector;
+    typedef typename Dune::BitSetVector<1> BitVector;
+
+    typedef typename GridType::LeafGridView GridView;
+    typedef typename Dune::FieldVector<double, dim> DomainType;
+    typedef typename Dune::FieldVector<double, 1> RangeType;
+
+
+    const GridView gridView = grid.leafView();
+
+    Matrix A;
+    constructPQ1Pattern(gridView, A);
+    A=0.0;
+    assemblePQ1Stiffness(gridView, A);
+
+    Vector rhs(A.N());
+    rhs = 0;
+    ConstantFunction<DomainType, RangeType> f(50);
+    assemblePQ1RHS(gridView, rhs, f);
+
+    Vector u(A.N());
+    u = 0;
+
+    BitVector ignore(A.N());
+    ignore.unsetAll();
+    markBoundaryDOFs(gridView, ignore);
+
+
+    solveObstacleProblemByMMGSolver(grid, A, u, rhs, ignore);
+
+    if (fileName!="")
+    {
+        typename Dune::VTKWriter<GridView> vtkWriter(gridView);
+        vtkWriter.addVertexData(u, "solution");
+        vtkWriter.write(fileName);
+    }
+
+    return passed;
+}
+
+
+template <int dim>
+bool checkWithYaspGrid(int refine, const std::string fileName="")
+{
+    bool passed = true;
+
+    typedef Dune::YaspGrid<dim> GridType;
+
+    typename Dune::FieldVector<typename GridType::ctype,dim> L(1.0);
+    typename Dune::FieldVector<int,dim> s(1);
+    typename Dune::FieldVector<bool,dim> periodic(false);
+
+    GridType grid(L, s, periodic, 0);
+
+    for (int i = 0; i < refine; ++i)
+        grid.globalRefine(1);
+
+    std::cout << "Testing with YaspGrid<" << dim << ">" << std::endl;
+    std::cout << "Number of levels: " << (grid.maxLevel()+1) << std::endl;
+    std::cout << "Number of leaf nodes: " << grid.leafView().size(dim) << std::endl;
+
+    passed = passed and checkWithGrid(grid, fileName);
+
+
+    return passed;
+}
+
+
+
+
+int main(int argc, char** argv) try
+{
+    bool passed(true);
+
+
+    int refine1d = 16;
+    int refine2d = 8;
+    int refine3d = 5;
+
+    if (argc>1)
+        std::istringstream(argv[1]) >> refine1d;
+    if (argc>2)
+        std::istringstream(argv[2]) >> refine2d;
+    if (argc>3)
+        std::istringstream(argv[3]) >> refine3d;
+
+    passed = passed and checkWithYaspGrid<1>(refine1d, "obstacletnnmgtest_yasp_1d_solution");
+    passed = passed and checkWithYaspGrid<2>(refine2d, "obstacletnnmgtest_yasp_2d_solution");
+    passed = passed and checkWithYaspGrid<3>(refine3d, "obstacletnnmgtest_yasp_3d_solution");
+
+    return passed ? 0 : 1;
+}
+catch (Dune::Exception e) {
+
+    std::cout << e << std::endl;
+
+}
+