diff --git a/dune/solvers/test/Makefile.am b/dune/solvers/test/Makefile.am index 1254b8f6c0f2c45536ff72a8cef3fbe90cafa33f..9175cce949728edec80fb661199102aee88c75ab 100644 --- a/dune/solvers/test/Makefile.am +++ b/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 diff --git a/dune/solvers/test/quadraticipoptsolvertest.cc b/dune/solvers/test/quadraticipoptsolvertest.cc new file mode 100644 index 0000000000000000000000000000000000000000..a7281954c56fccdbe8cdf145dafb8a16fa80a2f3 --- /dev/null +++ b/dune/solvers/test/quadraticipoptsolvertest.cc @@ -0,0 +1,196 @@ +#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; + +} +