Remove old example 'dgstokes'

test/navier-stokes/CMakeLists.txt

-add_subdirectory(howto)
+add_subdirectory(reference_program)
 
 dune_add_formcompiler_system_test(UFLFILE navierstokes_3d_dg_quadrilateral.ufl
                                   BASENAME navierstokes_3d_dg_quadrilateral

test/navier-stokes/howto/dgstokes.cc

-// -*- tab-width: 2; indent-tabs-mode: nil -*-
-/** \file
-    \brief Solve Stokes with DG method (stationary case).
-*/
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-#include <iostream>
-#include <vector>
-#include <map>
-#include <string>
-#include <random>
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
-#include <dune/istl/bvector.hh>
-#include <dune/istl/operators.hh>
-#include <dune/istl/solvers.hh>
-#include <dune/istl/preconditioners.hh>
-#include <dune/istl/io.hh>
-
-#include "dune/pdelab/finiteelementmap/qkdg.hh"
-#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
-#include<dune/pdelab/gridfunctionspace/subspace.hh>
-#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
-#include <dune/pdelab/gridfunctionspace/vtk.hh>
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-#include <dune/pdelab/gridfunctionspace/interpolate.hh>
-#include <dune/pdelab/localoperator/dgnavierstokes.hh>
-#include <dune/pdelab/backend/istl.hh>
-#include <dune/pdelab/finiteelementmap/monomfem.hh>
-#include <dune/pdelab/common/function.hh>
-#include <dune/pdelab/common/vtkexport.hh>
-
-#include "navierstokes_initial.hh"
-
-#define USE_SUPER_LU
-#define MAKE_VTK_OUTPUT
-
-//===============================================================
-// Problem setup and solution
-//===============================================================
-
-// generate a composite P2/P1 function and output it
-template<typename GV, typename RF, int vOrder, int pOrder>
-void stokes (const GV& gv, const Dune::ParameterTree& configuration, std::string filename)
-{
-  // <<<1>>> constants and types
-  using ES = Dune::PDELab::AllEntitySet<GV>;
-  ES es(gv);
-  typedef typename ES::Grid::ctype DF;
-  static const unsigned int dim = ES::dimension;
-  Dune::Timer watch;
-  std::cout << "=== Initialize" << std::endl;
-
-  // <<<2>>> Make grid function space
-  watch.reset();
-  // typedef Dune::PDELab::MonomLocalFiniteElementMap<DF,RF,dim,vOrder> vFEM;
-  // typedef Dune::PDELab::MonomLocalFiniteElementMap<DF,RF,dim,pOrder> pFEM;
-  // vFEM vFem(Dune::GeometryType(Dune::GeometryType::cube,dim));
-  // pFEM pFem(Dune::GeometryType(Dune::GeometryType::cube,dim));
-  typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, 2, 2> vFEM;
-  typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, 1, 2> pFEM;
-  vFEM vFem;
-  pFEM pFem;
-  // DOFs per cell
-  static const unsigned int vBlockSize = Dune::MonomImp::Size<dim,vOrder>::val;
-  static const unsigned int pBlockSize = Dune::MonomImp::Size<dim,pOrder>::val;
-
-
-  // typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,vBlockSize> VVectorBackend;
-  // typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,pBlockSize> PVectorBacken
-  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> VVectorBackend;
-  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> PVectorBackend;
-  typedef Dune::PDELab::istl::VectorBackend<> VelocityVectorBackend;
-
-
-#if 1
-  // this creates a flat backend (i.e. blocksize == 1)
-  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> VectorBackend;
-#else
-  // this creates a backend with static blocks matching the size of the LFS
-  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> VectorBackend;
-#endif
-  // velocity
-  Dune::dinfo << "--- v^dim" << std::endl;
-
-
-  // typedef Dune::PDELab::EntityBlockedOrderingTag VelocityOrderingTag;
-  typedef Dune::PDELab::LexicographicOrderingTag VelocityOrderingTag;
-
-
-  typedef Dune::PDELab::VectorGridFunctionSpace<
-    ES,vFEM,dim,
-    VelocityVectorBackend,
-    VVectorBackend,
-    Dune::PDELab::NoConstraints,
-    VelocityOrderingTag
-    > velocityGFS;
-  velocityGFS velocityGfs(es,vFem);
-  velocityGfs.name("v");
-  // p
-  Dune::dinfo << "--- p" << std::endl;
-  typedef Dune::PDELab::GridFunctionSpace<ES,pFEM,
-                                          Dune::PDELab::NoConstraints, PVectorBackend> pGFS;
-  pGFS pGfs(es,pFem);
-  pGfs.name("p");
-  // GFS
-  Dune::dinfo << "--- v^dim,p" << std::endl;
-
-
-  // typedef Dune::PDELab::EntityBlockedOrderingTag StokesOrderingTag;
-  typedef Dune::PDELab::LexicographicOrderingTag StokesOrderingTag;
-
-
-  typedef Dune::PDELab::CompositeGridFunctionSpace<
-    VectorBackend, StokesOrderingTag,
-    velocityGFS, pGFS> GFS;
-  GFS gfs(velocityGfs, pGfs);
-  using namespace Dune::Indices;
-  velocityGfs.child(_0).name("velocity_0");
-  velocityGfs.child(_1).name("velocity_1");
-  pGfs.name("pressure");
-  gfs.name("test");
-  gfs.update();
-  std::cout << "gfs with " << gfs.size() << " dofs generated  "<< std::endl;
-  std::cout << "=== function space setup " <<  watch.elapsed() << " s" << std::endl;
-
-  // <<<3>>> Make coefficient Vector and initialize it from a function
-  using V = Dune::PDELab::Backend::Vector<GFS,RF>;
-  V x(gfs);
-
-  typedef ZeroVectorFunction<ES,RF,dim> FType;
-  FType f(es);
-  typedef BCTypeParamGlobalDirichlet BType;
-  BType b;
-  typedef HagenPoiseuilleVelocityBox<ES,RF,dim> VType;
-  VType v(es);
-  typedef ZeroScalarFunction<ES,RF> PType;
-  PType p(es);
-
-  // <<<4>>> Make grid Function operator
-  watch.reset();
-  typedef Dune::PDELab::DefaultInteriorPenalty<RF> PenaltyTerm;
-
-  typedef Dune::PDELab::DGNavierStokesParameters<ES,RF,FType,BType,VType,PType,true,false,PenaltyTerm>
-    LocalDGOperatorParameters;
-  LocalDGOperatorParameters lop_params(configuration.sub("parameters"),f,b,v,p);
-  typedef Dune::PDELab::DGNavierStokes<LocalDGOperatorParameters> LocalDGOperator;
-  const int superintegration_order = 0;
-  LocalDGOperator lop(lop_params,superintegration_order);
-
-  typedef Dune::PDELab::EmptyTransformation C;
-  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
-  MBE mbe(75); // Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
-  typedef Dune::PDELab::GridOperator
-    <GFS,GFS,LocalDGOperator,MBE,RF,RF,RF,C,C> GOS;
-  GOS gos(gfs,gfs,lop,mbe);
-
-  std::cout << "=== grid operator space setup " <<  watch.elapsed() << " s" << std::endl;
-
-  typedef typename GOS::Jacobian M;
-  watch.reset();
-  M m(gos);
-  std::cout << m.patternStatistics() << std::endl;
-  std::cout << "=== matrix setup " <<  watch.elapsed() << " s" << std::endl;
-  m = 0.0;
-  watch.reset();
-  gos.jacobian(x,m);
-  std::cout << "=== jacobian assembly " <<  watch.elapsed() << " s" << std::endl;
-
-  // std::ofstream matrix("Matrix");
-  // Dune::printmatrix(matrix, m.base(), "M", "r", 6, 3);
-
-  // evaluate residual w.r.t initial guess
-  V r(gfs);
-  r = 0.0;
-  watch.reset();
-  x = 0.0;
-  gos.residual(x,r);
-  std::cout << "=== residual evaluation " <<  watch.elapsed() << " s" << std::endl;
-
-  bool verbose = true;
-
-  using Dune::PDELab::Backend::native;
-  typedef Dune::PDELab::Backend::Native<M> ISTLM;
-  typedef Dune::PDELab::Backend::Native<V> ISTLV;
-#ifdef USE_SUPER_LU // use lu decomposition as solver
-#if HAVE_SUPERLU
-  // make ISTL solver
-  Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(native(m));
-  Dune::SuperLU<ISTLM> solver(native(m), verbose?1:0);
-  Dune::InverseOperatorResult stat;
-#else
-#error No superLU support, please install and configure it.
-#endif
-#else // Use iterative solver
-  // make ISTL solver
-  Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(native(m));
-  Dune::SeqILU0<ISTLM,ISTLV,ISTLV> ilu0(native(m),1.0);
-  Dune::BiCGSTABSolver<ISTLV> solver(opa,ilu0,1E-10,20000, verbose?2:1);
-  Dune::InverseOperatorResult stat;
-#endif
-
-  // solve the jacobian system
-  r *= -1.0; // need -residual
-  x = r;
-  solver.apply(native(x),native(r),stat);
-
-
-
-  if (false) {
-    using Dune::PDELab::Backend::native;
-    V r(x);
-    // Setup random input
-    std::size_t seed = 0;
-    auto rng = std::mt19937_64(seed);
-    auto dist = std::uniform_real_distribution<>(-1., 1.);
-    for (auto& v : x)
-      v = dist(rng);
-    r=0.0;
-    gos.residual(x, r);
-    Dune::printvector(std::cout, native(r), "residual vector", "row");
-  }
-  if (true) {
-    // Setup random input
-    std::size_t seed = 0;
-    auto rng = std::mt19937_64(seed);
-    auto dist = std::uniform_real_distribution<>(1., 10.);
-    for (auto& v : x)
-      v = dist(rng);
-    //using M = typename GO::Traits::Jacobian;
-    M m(gos);
-    gos.jacobian(x, m);
-    using Dune::PDELab::Backend::native;
-    Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
-  }
-
-
-
-  //    #ifdef MAKE_VTK_OUTPUT
-  // output grid function with SubsamplingVTKWriter
-  Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,2);
-  Dune::PDELab::addSolutionToVTKWriter(vtkwriter, gfs, x);
-  vtkwriter.write(filename,Dune::VTK::appendedraw);
-  //    #endif
-}
-
-//===============================================================
-// Main program with grid setup
-//===============================================================
-
-int main(int argc, char** argv)
-{
-  try{
-    //Maybe initialize Mpi
-    Dune::MPIHelper::instance(argc, argv);
-
-    int x=2;
-    int y=2;
-    int z=2;
-    if (argc > 1)
-      x = atoi(argv[1]);
-    if (argc > 2)
-      y = atoi(argv[2]);
-    if (argc > 3)
-      z = atoi(argv[3]);
-
-    Dune::ParameterTree configuration;
-    const std::string config_filename("dgstokes.ini");
-    std::cout << "Reading ini-file \""<< config_filename
-              << "\"" << std::endl;
-
-    Dune::ParameterTreeParser::readINITree(config_filename, configuration);
-
-    // YaspGrid P2/P1 2D test
-    if(1){
-      // make grid
-      const int dim = 2;
-      Dune::FieldVector<double,dim> L(1.0);
-      Dune::array<int,dim> N;
-      N[0] = x; N[1] = y;
-      std::bitset<dim> B(false);
-      Dune::YaspGrid<dim> grid(L,N,B,0);
-
-      // get view
-      typedef Dune::YaspGrid<dim>::LeafGridView GV;
-      const GV gv=grid.leafGridView();
-
-      // solve problem
-      Dune::dinfo.push(false);
-      stokes<GV,double,2,1>(gv,configuration,"dgstokes-2D-2-1");
-    }
-
-
-    // YaspGrid P2/P1 3D test
-#if 0
-    {
-      // make grid
-      const int dim = 3;
-      Dune::FieldVector<double,dim> L(1.0);
-      Dune::array<int,dim> N;
-      N[0] = x; N[1] = y; N[2] = z;
-      std::bitset<dim> B(false);
-      Dune::YaspGrid<dim> grid(L,N,B,0);
-
-      // get view
-      typedef Dune::YaspGrid<dim>::LeafGridView GV;
-      const GV gv=grid.leafGridView();
-
-      // solve problem
-      stokes<GV,double,2,1>(gv,configuration,"dgstokes-3D-2-1");
-    }
-#endif
-
-    // test passed
-    return 0;
-  }
-  catch (Dune::Exception &e){
-    std::cerr << "Dune reported error: " << e << std::endl;
-    return 1;
-  }
-  catch (...){
-    std::cerr << "Unknown exception thrown!" << std::endl;
-    return 1;
-  }
-}

test/navier-stokes/howto/dgstokes.ini

-# Configuration file for the example programs - dgstokes, dgnavierstokesvelvecfem -
-
-[parameters]
-rho = 1.0
-mu = 1.0
-[parameters.dg]
-epsilon = -1
-sigma = 1.0
-beta = 1.0

test/navier-stokes/howto/navierstokes_initial.hh

-#ifndef CG_STOKES_INITIAL_HH
-#define CG_STOKES_INITIAL_HH
-
-//===============================================================
-// Define parameter functions f,g,j and \partial\Omega_D/N
-//===============================================================
-
-/** \brief Global Dirichlet boundary conditions
-
- */
-
-// constraints parameter class for selecting boundary condition type
-class BCTypeParamGlobalDirichlet
-{
-public :
-  typedef Dune::PDELab::StokesBoundaryCondition BC;
-
-  struct Traits
-  {
-    typedef BC::Type RangeType;
-  };
-
-  BCTypeParamGlobalDirichlet() {}
-
-  template<typename I>
-  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
-                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
-                        BC::Type& y) const
-  {
-    y = BC::VelocityDirichlet;
-  }
-
-};
-
-/** \brief Boundary conditions for Hagen-Poiseuille flow.
-
- */
-
-// constraints parameter class for selecting boundary condition type
-class BCTypeParamHagenPoiseuille
-{
-public:
-  typedef Dune::PDELab::StokesBoundaryCondition BC;
-
-  struct Traits
-  {
-    typedef BC::Type RangeType;
-  };
-
-  BCTypeParamHagenPoiseuille() {}
-
-  template<typename I>
-  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
-                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
-                        BC::Type& y) const
-  {
-    Dune::FieldVector<typename I::ctype, I::dimension>
-      xg = intersection.geometry().global( coord );
-    if( xg[0] > 1.0-1e-6 )
-      y = BC::DoNothing;
-    else
-      y = BC::VelocityDirichlet;
-  }
-};
-
-/** \brief Boundary conditions for two dimensional flow around a cylinder obstacle.
-
- * The settings come from the two dimensional Turek benchmark channel.
- * Compare with the Article: "Benchmark Computations of Laminar Flow Around a Cylinder"
- * by M. Schaefer and S. Turek.
-
- */
-
-// constraints parameter class for selecting boundary condition type
-class BCTypeParamTU
-{
-public:
-  typedef Dune::PDELab::StokesBoundaryCondition BC;
-
-  struct Traits
-  {
-    typedef BC::Type RangeType;
-  };
-
-  BCTypeParamTU() {}
-
-  template<typename I>
-  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
-                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
-                        BC::Type& y) const
-  {
-    Dune::FieldVector<typename I::ctype, I::dimension>
-      xg = intersection.geometry().global( coord );
-    if( xg[0] > 2.2-1e-6 )
-      y = BC::DoNothing;
-    else
-      y = BC::VelocityDirichlet;
-  }
-
-  template<typename Time>
-  void setTime(Time)
-  {}
-};
-
-/** \brief Boundary conditions for three dimensional flow around a cylinder obstacle
-
- * The settings come from the three dimensional Turek benchmark channel.
- * Compare with the Article: "Benchmark Computations of 3D Laminar Flow Around a Cylinder with CFX, OpenFOAM and FEATFLOW"
- * by E. Bayraktar, O. Mierka and S. Turek.
-
- */
-
-// constraints parameter class for selecting boundary condition type
-class BCTypeParamTU3D
-{
-public:
-  typedef Dune::PDELab::StokesBoundaryCondition BC;
-
-  struct Traits
-  {
-    typedef BC::Type RangeType;
-  };
-
-  BCTypeParamTU3D() {}
-
-  template<typename I>
-  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
-                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
-                        BC::Type& y) const
-  {
-    Dune::FieldVector<typename I::ctype, I::dimension>
-      xg = intersection.geometry().global( coord );
-    if( xg[0] > 2.5-1e-6 )
-      y = BC::DoNothing;
-    else
-      y = BC::VelocityDirichlet;
-  }
-};
-
-/** \brief Boundary conditions for Backward-Facing Step.
-
- */
-
-// constraints parameter class for selecting boundary condition type
-class BCTypeParamLU
-{
-public:
-  typedef Dune::PDELab::StokesBoundaryCondition BC;
-
-  struct Traits
-  {
-    typedef BC::Type RangeType;
-  };
-
-  BCTypeParamLU() {}
-
-  template<typename I>
-  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
-                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
-                        BC::Type& y) const
-  {
-    Dune::FieldVector<typename I::ctype, I::dimension>
-      xg = intersection.geometry().global( coord );
-    if( xg[0] > 5.0-1e-6 )
-      y = BC::DoNothing;
-    else
-      y = BC::VelocityDirichlet;
-  }
-
-  template<typename Time>
-  void setTime(Time)
-  {}
-};
-
-/** \brief Dirichlet velocity function for Hagen-Poiseuille flow in a box.
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-    \tparam dim Dimension of the problem.
-
-    * The Hagen-Poiseuille flow is solved on the dim-dimensional unitcube.
-    * The maximum inflow velocity is normalized to 1.
-    * In three dimensions, the two dimensional flow is trivially extended
-    * to the z-direction.
-
-    */
-
-template<typename GV, typename RF, int dim>
-class HagenPoiseuilleVelocityBox :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
-  HagenPoiseuilleVelocityBox<GV,RF,dim> >
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, HagenPoiseuilleVelocityBox<GV,RF,dim> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  HagenPoiseuilleVelocityBox(const GV & gv) : BaseT(gv) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0;
-    y[0] = x[1]*(1.0-x[1]);
-    y[0] *= 4.0;
-  }
-
-};
-
-/** \brief Dirichlet velocity function for the flow around a cylinder obstacle.
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-    \tparam dim The dimension of the problem.
-    \param[in] meanflow Maximum inflow velocity.
-
-    * This function is used to model the inflow.
-    * It works both in two and three dimensions.
-
-    */
-
-template<typename GV, typename RF, int dim>
-class TUVelocity :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
-  TUVelocity<GV,RF,dim> >
-{
-  RF time;
-  const RF meanflow_;
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TUVelocity<GV,RF,dim> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  TUVelocity(const GV & gv, const RF& meanflow) :
-    BaseT(gv)
-    , meanflow_(meanflow)
-  {
-    time = 0.0;
-  }
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0.0;
-    RF r = 1.0;
-
-    for(int i=1; i<dim; i++)
-      r *= 4.0*x[i]*(0.41-x[i])/(0.41*0.41);
-
-    if(x[0] < 1e-6) {
-      y[0] = r*meanflow_;
-      if(time <= 4.0)
-        y[0] *= sin(M_PI*time/8.0);
-    }
-  }
-
-  template <typename T>
-  void setTime(T t){
-    time = t;
-  }
-
-};
-
-/** \brief Dirichlet velocity function for the Backward-Facing Step.
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-    \tparam dim The dimension of the problem.
-
-    * The inflow velocity is normalized to 1.
-
-    */
-
-template<typename GV, typename RF, int dim>
-class LUVelocity :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
-  LUVelocity<GV,RF,dim> >
-{
-private:
-  RF time;
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,LUVelocity<GV,RF,dim> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  LUVelocity(const GV & gv) : BaseT(gv) {
-    time = 0.0;
-  }
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    RF r = 0;
-
-    y = 0.0;
-    r += x[1]*(1.0-x[1]);
-
-    if(x[0] < -1.0+1e-6) {
-      y[0] = 4*r;
-      if(time <= 2.0)
-        y[0] *= sin(M_PI/4*time);
-    }
-  }
-
-  template <typename T>
-  void setTime(T t){
-    time = t;
-  }
-
-};
-
-/** \brief Dirichlet velocity function for three dimensional Hagen-Poiseuille flow.
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-
-    * The Hagen-Poiseuille flow is solved in a pipe pointing in x-direction.
-    * The cross section in the yz-plane is a circle with midpoint at the origin
-    * and with radius 0.5.
-    * The inflow velocity is normalized to 1.
-
-    */
-
-template<typename GV, typename RF>
-class HagenPoiseuilleVelocityCylindrical :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,3>,
-  HagenPoiseuilleVelocityCylindrical<GV,RF> >
-{
-public:
-  enum {dim = 3};
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, HagenPoiseuilleVelocityCylindrical<GV,RF> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  HagenPoiseuilleVelocityCylindrical(const GV & gv) : BaseT(gv)
-  {
-    assert(GV::dimension == dim);
-  }
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-
-    RF r = 0;
-    for(int i=1; i<dim; ++i){
-      r += (x[i])*(x[i]);
-      y[i] = 0;
-    }
-    r = sqrt(r);
-    y[0] = 0.25 - r*r;
-    y[0] *= 4;
-
-    if(y[0]<0)
-      y[0] = 0;
-  }
-
-};
-
-
-template<typename GV, typename RF, std::size_t dim_range>
-class ZeroVectorFunction :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,
-  ZeroVectorFunction<GV,RF,dim_range> >,
-  public Dune::PDELab::InstationaryFunctionDefaults
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroVectorFunction> BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  ZeroVectorFunction(const GV & gv) : BaseT(gv) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y=0;
-  }
-};
-
-template<typename GV, typename RF>
-class ZeroScalarFunction
-  : public ZeroVectorFunction<GV,RF,1>
-{
-public:
-
-  ZeroScalarFunction(const GV & gv) : ZeroVectorFunction<GV,RF,1>(gv) {}
-
-};
-
-/** \brief Do-nothing boundary function for the Hagen-Poiseuille flow
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-
-*/
-
-// function for defining the flux boundary condition
-template<typename GV, typename RF>
-class HagenPoiseuilleZeroFlux
-  : public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
-  HagenPoiseuilleZeroFlux<GV,RF> >
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,HagenPoiseuilleZeroFlux<GV,RF> > BaseT;
-
-  HagenPoiseuilleZeroFlux (const GV& gv) : BaseT(gv) {}
-
-  inline void evaluateGlobal (const typename Traits::DomainType& x,
-                              typename Traits::RangeType& y) const
-  {
-    y = 0;
-  }
-};
-
-/** \brief Do-nothing boundary function for the Backward-Facing Step and the flow around a cylinder
-
-    \tparam GV The underlying grid view.
-    \tparam RF The type of the range field.
-    \param[in] p_ The pressure at the inflow.
-    \param[in] l_ The length of the tube.
-    \param[in] o_ The coordinate referring to the beginning of the tube.
-    \param[in] d_ The canonical direction the tube is oriented.
-
-*/
-
-// function for defining the flux boundary condition
-template<typename GV, typename RF>
-class PressureDropFlux
-  : public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
-  PressureDropFlux<GV,RF> >
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,PressureDropFlux<GV,RF> > BaseT;
-
-private:
-  typedef typename Traits::DomainFieldType DFT;
-
-  const DFT pressure;
-  const DFT length;
-  const DFT origin;
-  const int direction;
-
-public:
-  PressureDropFlux (const GV& gv, const RF p_, const RF l_, const RF o_, const int d_)
-    : BaseT(gv), pressure(p_), length(l_), origin(o_), direction(d_)
-  {
-#ifndef NDEBUG
-    const int dim = GV::dimension;
-    assert(direction >=0 && direction <dim);
-#endif
-  }
-
-  inline void evaluateGlobal (const typename Traits::DomainType& x,
-                              typename Traits::RangeType& y) const
-  {
-    if(x[direction]-origin < 1e-6)
-      y = pressure;
-    else if(x[direction]-origin > length-1e-6)
-      y = 0.;
-  }
-
-  template<typename Time>
-  void setTime(Time)
-  {}
-};
-
-
-#endif

test/navier-stokes/howto/CMakeLists.txt → test/navier-stokes/reference_program/CMakeLists.txt

-add_executable(dgstokes dgstokes.cc)
 add_executable(taylor-green taylor-green.cc)
 
-dune_symlink_to_source_files(FILES dgstokes.ini taylor-green.ini)
+dune_symlink_to_source_files(FILES taylor-green.ini)