add a rudimentary test for the QuadraticIPOptSolver class

[[Imported from SVN: r7294]]

dune/solvers/test/Makefile.am

@@ -3,6 +3,7 @@
 TESTS = genericvectortoolstest \
         lowrankoperatortest \
         nulloperatortest \
+        quadraticipoptsolvertest \
         sumoperatortest \
         obstacletnnmgtest
 
@@ -37,14 +38,22 @@ nulloperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
 nulloperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
 nulloperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
 
+quadraticipoptsolvertest_SOURCES = quadraticipoptsolvertest.cc
+# quadraticipoptsolvertest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) $(IPOPT_CPPFLAGS)
+# quadraticipoptsolvertest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) $(IPOPT_LDFLAGS)
+# quadraticipoptsolvertest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) $(IPOPT_LIBS)
+quadraticipoptsolvertest_CPPFLAGS = $(COMMON_CPPFLAGS) $(IPOPT_CPPFLAGS)
+quadraticipoptsolvertest_LDADD = $(COMMON_LDADD) $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
+quadraticipoptsolvertest_LDFLAGS = $(COMMON_LDFLAGS) $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
+
 sumoperatortest_SOURCES = sumoperatortest.cc
 sumoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
 sumoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
 sumoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
 
 obstacletnnmgtest_SOURCES = obstacletnnmgtest.cc
-obstacletnnmgtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) #${IPOPT_CPPFLAGS}
-obstacletnnmgtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) #${IPOPT_LIBS}
-obstacletnnmgtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) #${IPOPT_LDFLAGS}
+obstacletnnmgtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) ${IPOPT_CPPFLAGS}
+obstacletnnmgtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) ${IPOPT_LIBS}
+obstacletnnmgtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) ${IPOPT_LDFLAGS}
 
 include $(top_srcdir)/am/global-rules

dune/solvers/test/quadraticipoptsolvertest.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/solvers/quadraticipopt.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_;
+};
+
+
+#if 0
+template<class GridType, class MatrixType, class VectorType>
+void solveObstacleProblemByTNNMG(const GridType& grid, const MatrixType& mat, VectorType& x, const VectorType& rhs, int maxIterations=100, double tolerance=1.0e-10)
+{
+    // bit vector marking dofs to ignore
+    typedef typename ObstacleTNNMGStep<MatrixType, VectorType>::BitVector BitVector;
+    BitVector ignore(rhs.size());
+    ignore.unsetAll();
+
+    solveObstacleProblemByTNNMG(grid, mat, x, rhs, ignore, maxIterations, tolerance);
+}
+#endif
+
+template<class GridType, class MatrixType, class VectorType, class BitVector>
+void solveObstacleProblemByQuadraticIPOptSolver(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;
+
+    const int blockSize = VectorType::block_type::dimension;
+
+    // create double obstacle constraints
+#if 0
+    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 solver
+    Solver solver;
+
+    // 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);
+
+
+    solveObstacleProblemByQuadraticIPOptSolver(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;
+
+}
+