<commit_before>﻿// This file is part of the dune-gdt project:
//   http://users.dune-project.org/projects/dune-gdt
// Copyright holders: <NAME>
// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

#ifndef DUNE_GDT_TEST_ELLIPTIC_CG_DISCRETIZATION_HH
#define DUNE_GDT_TEST_ELLIPTIC_CG_DISCRETIZATION_HH

#if !HAVE_DUNE_PDELAB
# error "This one requires dune-pdelab!"
#endif

#include <memory>
#include <vector>
#include <type_traits>

#include <dune/common/timer.hh>

#include <dune/stuff/common/memory.hh>
#include <dune/stuff/grid/boundaryinfo.hh>
#include <dune/stuff/la/container/eigen.hh>
#include <dune/stuff/la/solver.hh>
#include <dune/stuff/functions/interfaces.hh>
#include <dune/stuff/common/convergence-study.hh>
#include <dune/stuff/functions/combined.hh>

#include <dune/gdt/spaces/continuouslagrange/pdelab.hh>
#include <dune/gdt/operators/elliptic.hh>
#include <dune/gdt/functionals/l2.hh>
#include <dune/gdt/assembler/system.hh>
#include <dune/gdt/spaces/constraints.hh>
#include <dune/gdt/products/l2.hh>
#include <dune/gdt/products/h1.hh>
#include <dune/gdt/operators/prolongations.hh>
#include <dune/gdt/operators/projections.hh>

#include "elliptic-testcases.hh"


namespace EllipticCG {


template< class GridViewType,
          int polynomialOrder,
          class MatrixImp = Dune::Stuff::LA::EigenRowMajorSparseMatrix< double >,
          class VectorImp = Dune::Stuff::LA::EigenDenseVector< double > >
class Discretization
{
public:
  static const unsigned int             dimDomain = GridViewType::dimension;
  typedef typename GridViewType::ctype  DomainFieldType;

  static const unsigned int dimRange = 1;
  typedef double            RangeFieldType;

  static const unsigned int polOrder = polynomialOrder;

  typedef Dune::Stuff::Grid::BoundaryInfoInterface< typename GridViewType::Intersection > BoundaryInfoType;
  typedef Dune::Stuff::LocalizableFunctionInterface
      < typename GridViewType::template Codim< 0 >::Entity, DomainFieldType, dimDomain, RangeFieldType, dimRange >
    FunctionType;

  typedef MatrixImp MatrixType;
  typedef VectorImp VectorType;

  typedef Dune::GDT::Spaces::ContinuousLagrange::PdelabBased
      < GridViewType, polOrder, RangeFieldType, dimRange > SpaceType;

  typedef Dune::GDT::DiscreteFunction< SpaceType, VectorType >      DiscreteFunctionType;
  typedef Dune::GDT::ConstDiscreteFunction< SpaceType, VectorType > ConstDiscreteFunctionType;

  Discretization(const std::shared_ptr< const GridViewType >& gp,
                 const BoundaryInfoType& info,
                 const FunctionType& diff,
                 const FunctionType& forc,
                 const FunctionType& dir,
                 const FunctionType& neu)
    : space_(gp)
    , boundary_info_(info)
    , diffusion_(diff)
    , force_(forc)
    , dirichlet_(dir)
    , neumann_(neu)
    , is_assembled_(false)
    , system_matrix_(0, 0)
    , rhs_vector_(0)
    , dirichlet_shift_vector_(0)
  {}

  const SpaceType& space() const
  {
    return space_;
  }

  VectorType create_vector() const
  {
    return VectorType(space_.mapper().size());
  }

  void assemble() const
  {
    using namespace Dune;
    using namespace Dune::GDT;
    if (!is_assembled_) {
      // create the containers (use the sparsity pattern of the operator)
      typedef Operators::EllipticCG< FunctionType, MatrixType, SpaceType > EllipticOperatorType;
      system_matrix_ = MatrixType(space_.mapper().size(),
                                  space_.mapper().size(),
                                  EllipticOperatorType::pattern(space_));
      rhs_vector_ = VectorType(space_.mapper().size());
      dirichlet_shift_vector_ = VectorType(space_.mapper().size());

      // define the lhs operator and the rhs functionals
      EllipticOperatorType elliptic_operator(diffusion_, system_matrix_, space_);
      typedef GDT::Functionals::L2Volume< FunctionType, VectorType, SpaceType > L2VolumeFunctionalType;
      L2VolumeFunctionalType force_functional(force_, rhs_vector_, space_);
      typedef GDT::Functionals::L2Face< FunctionType, VectorType, SpaceType > L2FaceFunctionalType;
      L2FaceFunctionalType neumann_functional(neumann_, rhs_vector_, space_);
      // project the dirichlet boundary values
      DiscreteFunctionType dirichlet_projection(space_, dirichlet_shift_vector_);
      typedef Operators::DirichletProjectionLocalizable< GridViewType, FunctionType, DiscreteFunctionType >
          DirichletProjectionOperator;
      DirichletProjectionOperator dirichlet_projection_operator(*(space_.grid_view()),
                                                                boundary_info_,
                                                                dirichlet_,
                                                                dirichlet_projection);
      // assemble everything
      SystemAssembler< SpaceType > grid_walker(space_);
      grid_walker.add(elliptic_operator);
      grid_walker.add(force_functional);
      grid_walker.add(neumann_functional, new ApplyOn::NeumannIntersections< GridViewType >(boundary_info_));
      grid_walker.add(dirichlet_projection_operator, new ApplyOn::BoundaryEntities< GridViewType >());
      grid_walker.walk();
      // substract the operators action on the dirichlet values
      auto tmp = rhs_vector_.copy();
      system_matrix_.mv(dirichlet_shift_vector_, tmp);
      rhs_vector_ -= tmp;
      // apply the dirichlet constraints
      Constraints::Dirichlet < typename GridViewType::Intersection, RangeFieldType >
        dirichlet_constraints(boundary_info_, space_.mapper().maxNumDofs(), space_.mapper().maxNumDofs());
      grid_walker.add(dirichlet_constraints, system_matrix_);
      grid_walker.add(dirichlet_constraints, rhs_vector_);
      grid_walker.walk();

      is_assembled_ = true;
    }
  } // ... assemble()

  bool assembled() const
  {
    return is_assembled_;
  }

  const MatrixType& system_matrix() const
  {
    return system_matrix_;
  }

  const VectorType& rhs_vector() const
  {
    return rhs_vector_;<commit_msg>[tests.elliptic-cg...] whitespace<commit_after>﻿// This file is part of the dune-gdt project:
//   http://users.dune-project.org/projects/dune-gdt
// Copyright holders: <NAME>
// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

#ifndef DUNE_GDT_TEST_ELLIPTIC_CG_DISCRETIZATION_HH
#define DUNE_GDT_TEST_ELLIPTIC_CG_DISCRETIZATION_HH

#if !HAVE_DUNE_PDELAB
# error "This one requires dune-pdelab!"
#endif

#include <memory>
#include <vector>
#include <type_traits>

#include <dune/common/timer.hh>

#include <dune/stuff/common/memory.hh>
#include <dune/stuff/grid/boundaryinfo.hh>
#include <dune/stuff/la/container/eigen.hh>
#include <dune/stuff/la/solver.hh>
#include <dune/stuff/functions/interfaces.hh>
#include <dune/stuff/common/convergence-study.hh>
#include <dune/stuff/functions/combined.hh>

#include <dune/gdt/spaces/continuouslagrange/pdelab.hh>
#include <dune/gdt/operators/elliptic.hh>
#include <dune/gdt/functionals/l2.hh>
#include <dune/gdt/assembler/system.hh>
#include <dune/gdt/spaces/constraints.hh>
#include <dune/gdt/products/l2.hh>
#include <dune/gdt/products/h1.hh>
#include <dune/gdt/operators/prolongations.hh>
#include <dune/gdt/operators/projections.hh>

#include "elliptic-testcases.hh"


namespace EllipticCG {


template< class GridViewType,
          int polynomialOrder,
          class MatrixImp = Dune::Stuff::LA::EigenRowMajorSparseMatrix< double >,
          class VectorImp = Dune::Stuff::LA::EigenDenseVector< double > >
class Discretization
{
public:
  static const unsigned int             dimDomain = GridViewType::dimension;
  typedef typename GridViewType::ctype  DomainFieldType;

  static const unsigned int dimRange = 1;
  typedef double            RangeFieldType;

  static const unsigned int polOrder = polynomialOrder;

  typedef Dune::Stuff::Grid::BoundaryInfoInterface< typename GridViewType::Intersection > BoundaryInfoType;
  typedef Dune::Stuff::LocalizableFunctionInterface
      < typename GridViewType::template Codim< 0 >::Entity, DomainFieldType, dimDomain, RangeFieldType, dimRange >
    FunctionType;

  typedef MatrixImp MatrixType;
  typedef VectorImp VectorType;

  typedef Dune::GDT::Spaces::ContinuousLagrange::PdelabBased
      < GridViewType, polOrder, RangeFieldType, dimRange > SpaceType;

  typedef Dune::GDT::DiscreteFunction< SpaceType, VectorType >      DiscreteFunctionType;
  typedef Dune::GDT::ConstDiscreteFunction< SpaceType, VectorType > ConstDiscreteFunctionType;

  Discretization(const std::shared_ptr< const GridViewType >& gp,
                 const BoundaryInfoType& info,
                 const FunctionType& diff,
                 const FunctionType& forc,
                 const FunctionType& dir,
                 const FunctionType& neu)
    : space_(gp)
    , boundary_info_(info)
    , diffusion_(diff)
    , force_(forc)
    , dirichlet_(dir)
    , neumann_(neu)
    , is_assembled_(false)
    , system_matrix_(0, 0)
    , rhs_vector_(0)
    , dirichlet_shift_vector_(0)
  {}

  const SpaceType& space() const
  {
    return space_;
  }

  VectorType create_vector() const
  {
    return VectorType(space_.mapper().size());
  }

  void assemble() const
  {
    using namespace Dune;
    using namespace Dune::GDT;
    if (!is_assembled_) {
      // create the containers (use the sparsity pattern of the operator)
      typedef Operators::EllipticCG< FunctionType, MatrixType, SpaceType > EllipticOperatorType;
      system_matrix_ = MatrixType(space_.mapper().size(),
                                  space_.mapper().size(),
                                  EllipticOperatorType::pattern(space_));
      rhs_vector_ = VectorType(space_.mapper().size());
      dirichlet_shift_vector_ = VectorType(space_.mapper().size());

      // define the lhs operator and the rhs functionals
      EllipticOperatorType elliptic_operator(diffusion_, system_matrix_, space_);
      typedef GDT::Functionals::L2Volume< FunctionType, VectorType, SpaceType > L2VolumeFunctionalType;
      L2VolumeFunctionalType force_functional(force_, rhs_vector_, space_);
      typedef GDT::Functionals::L2Face< FunctionType, VectorType, SpaceType > L2FaceFunctionalType;
      L2FaceFunctionalType neumann_functional(neumann_, rhs_vector_, space_);
      // project the dirichlet boundary values
      DiscreteFunctionType dirichlet_projection(space_, dirichlet_shift_vector_);
      typedef Operators::DirichletProjectionLocalizable< GridViewType, FunctionType, DiscreteFunctionType >
          DirichletProjectionOperator;
      DirichletProjectionOperator dirichlet_projection_operator(*(space_.grid_view()),
                                                                boundary_info_,
                                                                dirichlet_,
                                                                dirichlet_projection);
      // assemble everything
      SystemAssembler< SpaceType > grid_walker(space_);
      grid_walker.add(elliptic_operator);
      grid_walker.add(force_functional);
      grid_walker.add(neumann_functional, new ApplyOn::NeumannIntersections< GridViewType >(boundary_info_));
      grid_walker.add(dirichlet_projection_operator, new ApplyOn::BoundaryEntities< GridViewType >());
      grid_walker.walk();
      // substract the operators action on the dirichlet values
      auto tmp = rhs_vector_.copy();
      system_matrix_.mv(dirichlet_shift_vector_, tmp);
      rhs_vector_ -= tmp;
      // apply the dirichlet constraints
      Constraints::Dirichlet< typename GridViewType::Intersection, RangeFieldType >
        dirichlet_constraints(boundary_info_, space_.mapper().maxNumDofs(), space_.mapper().maxNumDofs());
      grid_walker.add(dirichlet_constraints, system_matrix_);
      grid_walker.add(dirichlet_constraints, rhs_vector_);
      grid_walker.walk();

      is_assembled_ = true;
    }
  } // ... assemble()

  bool assembled() const
  {
    return is_assembled_;
  }

  const MatrixType& system_matrix() const
  {
    return system_matrix_;
  }

  const VectorType& rhs_vector() const
  {
    return rhs_vector_;
{'fraction_non_alphanumeric': 0.09893730949831123, 'fraction_numerical': 0.002800889694373507, 'mean_word_length': 3.4698085419734905, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '956383', 'n_tokens_mistral': 3743, 'n_tokens_neox': 3549, 'n_words': 702}