cutfem_solver.hh 7.32 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
#ifndef DUNEURO_CUTFEM_SOLVER_HH
#define DUNEURO_CUTFEM_SOLVER_HH

#include <dune/udg/pdelab/cutfemmultiphaseoperator.hh>
#include <dune/udg/pdelab/multiphaseoperator.hh>
#include <dune/udg/pdelab/operator.hh>
#include <dune/udg/pdelab/subtriangulation.hh>

#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/stationary/linearproblem.hh>

#include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/convection_diffusion_udg_default_parameter.hh>
#include <duneuro/common/cutfem_gridoperator.hh>
#include <duneuro/common/cutfem_multi_phase_space.hh>
#include <duneuro/common/edge_norm_provider.hh>
Andreas Nüßing's avatar
Andreas Nüßing committed
17
#include <duneuro/common/kdtree.hh>
18
#include <duneuro/common/linear_problem_solver.hh>
19
#include <duneuro/common/penalty_flux_weighting.hh>
20
#include <duneuro/common/random.hh>
21
#include <duneuro/common/vector_initialization.hh>
22 23 24 25 26 27 28

namespace duneuro
{
  template <class ST, int comps, int degree, class P, class DF, class RF, class JF>
  struct CutFEMSolverTraits {
    using SubTriangulation = ST;
    using FundamentalGridView = typename ST::BaseT::GridView;
Andreas Nüßing's avatar
Andreas Nüßing committed
29 30
    using CoordinateFieldType = typename FundamentalGridView::ctype;
    using ElementSearch = KDTreeElementSearch<FundamentalGridView>;
31 32 33 34 35 36 37 38 39
    static const int dimension = FundamentalGridView::dimension;
    static const int compartments = comps;
    using Problem = P;
    using FunctionSpace = CutFEMMultiPhaseSpace<FundamentalGridView, RF, degree, compartments>;
    using DomainField = DF;
    using RangeField = RF;
    using DomainDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, DF>;
    using RangeDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, RF>;
    using EdgeNormProvider = MultiEdgeNormProvider;
40 41 42
    using PenaltyFluxWeighting = UnfittedDynamicPenaltyFluxWeights;
    using LocalOperator =
        ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
    using WrappedLocalOperator = Dune::UDG::CutFEMMultiPhaseLocalOperatorWrapper<LocalOperator>;
    // using WrappedLocalOperator = Dune::UDG::MultiPhaseLocalOperatorWrapper<LocalOperator>;
    using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<FundamentalGridView>;
    using MatrixBackend = Dune::PDELab::istl::BCRSMatrixBackend<>;
    using RawGridOperator =
        Dune::UDG::UDGGridOperator<typename FunctionSpace::GFS, typename FunctionSpace::GFS,
                                   WrappedLocalOperator, MatrixBackend, DF, RF, JF,
                                   UnfittedSubTriangulation>;
    using GridOperator = CutFEMGridOperator<RawGridOperator, ST, EdgeNormProvider>;
    using SolverBackend = Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GridOperator>;
    using LinearSolver = LinearProblemSolver<GridOperator, DomainDOFVector, RangeDOFVector>;
  };

  template <class ST, int compartments, int degree,
            class P = ConvectionDiffusion_UDG_DefaultParameter<typename ST::BaseT::GridView>,
            class DF = double, class RF = double, class JF = double>
  class CutFEMSolver
  {
  public:
    using Traits = CutFEMSolverTraits<ST, compartments, degree, P, DF, RF, JF>;

Andreas Nüßing's avatar
Andreas Nüßing committed
64 65
    CutFEMSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
                 std::shared_ptr<const typename Traits::ElementSearch> search,
66
                 const Dune::ParameterTree& config)
Andreas Nüßing's avatar
Andreas Nüßing committed
67 68 69
        : CutFEMSolver(subTriangulation, search,
                       std::make_shared<typename Traits::Problem>(
                           config.get<std::vector<double>>("conductivities")),
70 71 72 73
                       config)
    {
    }

Andreas Nüßing's avatar
Andreas Nüßing committed
74 75
    CutFEMSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
                 std::shared_ptr<const typename Traits::ElementSearch> search,
76 77 78
                 std::shared_ptr<typename Traits::Problem> problem,
                 const Dune::ParameterTree& config)
        : subTriangulation_(subTriangulation)
Andreas Nüßing's avatar
Andreas Nüßing committed
79
        , search_(search)
80 81 82
        , problem_(problem)
        , functionSpace_(subTriangulation_->gridView(), subTriangulation_)
        , edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
83
        , weighting_(config.get<std::string>("weights"))
84 85 86 87
        , localOperator_(
              *problem_, edgeNormProvider_, weighting_,
              ConvectionDiffusion_DG_Scheme::fromString(config.get<std::string>("scheme")),
              config.get<RF>("penalty"))
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
        , wrappedLocalOperator_(localOperator_)
        , unfittedSubTriangulation_(subTriangulation_->gridView(), *subTriangulation_)
        , rawGridOperator_(functionSpace_.getGFS(), functionSpace_.getGFS(),
                           unfittedSubTriangulation_, wrappedLocalOperator_,
                           typename Traits::MatrixBackend(2 * Traits::dimension + 1))
        , gridOperator_(rawGridOperator_, subTriangulation, edgeNormProvider_, config)
        , linearSolver_(gridOperator_, config)
    {
    }

    template <class SolverBackend>
    void solve(SolverBackend& solverBackend, const typename Traits::RangeDOFVector& rightHandSide,
               typename Traits::DomainDOFVector& solution, const Dune::ParameterTree& config,
               DataTree dataTree = DataTree())
    {
      Dune::Timer timer;
      randomize_uniform(Dune::PDELab::Backend::native(solution), DF(-1.0), DF(1.0));
      linearSolver_.apply(solverBackend, solution, rightHandSide, config, dataTree);
      dataTree.set("time", timer.elapsed());
    }

    template <class SolverBackend>
    void solve(SolverBackend& solverBackend, typename Traits::DomainDOFVector& solution,
               const Dune::ParameterTree& config, DataTree dataTree = DataTree())
    {
      Dune::Timer timer;
114 115 116
      initialize(Dune::PDELab::Backend::native(solution), config.hasSub("initialization") ?
                                                              config.sub("initialization") :
                                                              Dune::ParameterTree());
117 118 119 120 121 122 123 124 125
      linearSolver_.apply(solverBackend, solution, config, dataTree);
      dataTree.set("time", timer.elapsed());
    }

    const typename Traits::FunctionSpace& functionSpace() const
    {
      return functionSpace_;
    }

Andreas Nüßing's avatar
Andreas Nüßing committed
126
    std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation() const
127
    {
Andreas Nüßing's avatar
Andreas Nüßing committed
128
      return subTriangulation_;
129 130 131 132 133 134 135
    }

    typename Traits::Problem& problem()
    {
      return *problem_;
    }

Andreas Nüßing's avatar
Andreas Nüßing committed
136 137 138 139 140
    std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
    {
      return search_;
    }

141 142 143 144 145
    bool scaleToBBox() const
    {
      return false;
    }

146 147 148 149 150 151 152
#if HAVE_TBB
    tbb::mutex& functionSpaceMutex()
    {
      return fsMutex_;
    }
#endif

153
  private:
Andreas Nüßing's avatar
Andreas Nüßing committed
154 155
    std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation_;
    std::shared_ptr<const typename Traits::ElementSearch> search_;
156 157 158
    std::shared_ptr<typename Traits::Problem> problem_;
    typename Traits::FunctionSpace functionSpace_;
    typename Traits::EdgeNormProvider edgeNormProvider_;
159
    typename Traits::PenaltyFluxWeighting weighting_;
160 161 162 163 164 165
    typename Traits::LocalOperator localOperator_;
    typename Traits::WrappedLocalOperator wrappedLocalOperator_;
    typename Traits::UnfittedSubTriangulation unfittedSubTriangulation_;
    typename Traits::RawGridOperator rawGridOperator_;
    typename Traits::GridOperator gridOperator_;
    typename Traits::LinearSolver linearSolver_;
166 167 168 169

#if HAVE_TBB
    tbb::mutex fsMutex_;
#endif
170 171 172 173
  };
}

#endif // DUNEURO_CUTFEM_SOLVER_HH