adapt_solver.h

// This file is part of Hermes2D.
//
// Hermes2D is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// Hermes2D is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D.  If not, see <http://www.gnu.org/licenses/>.

#ifndef __H2D_ADAPT_SOLVER_H
#define __H2D_ADAPT_SOLVER_H

#include "adapt.h"
#include "error_calculator.h"
#include "../solver/linear_solver.h"
#include "../solver/newton_solver.h"
#include "../solver/picard_solver.h"
#include "../refinement_selectors/selector.h"
#include <unordered_set>

namespace Hermes
{
  namespace Hermes2D
  {
    class HERMES_API AdaptSolverCriterion
    {
    public:
      AdaptSolverCriterion();
      virtual bool done(double error, unsigned short iteration) = 0;
    };

    class HERMES_API AdaptSolverCriterionErrorThreshold : public AdaptSolverCriterion
    {
    public:
      AdaptSolverCriterionErrorThreshold(double error_threshold);
      virtual bool done(double error, unsigned short iteration);
      double error_threshold;
    };

    class HERMES_API AdaptSolverCriterionFixed : public AdaptSolverCriterion
    {
    public:
      AdaptSolverCriterionFixed(unsigned short refinement_levels);
      virtual bool done(double error, unsigned short iteration);
      unsigned short refinement_levels;
    };

    enum AdaptivityType
    {
      pAdaptivity,
      hAdaptivity,
      hpAdaptivity
    };

    template<typename Scalar, typename SolverType>
    class HERMES_API AdaptSolver :
      public Hermes::Mixins::TimeMeasurable,
      public Hermes::Mixins::Loggable,
      public Hermes::Mixins::StateQueryable
    {
    public:
      AdaptSolver(std::vector<SpaceSharedPtr<Scalar> > initial_spaces, WeakFormSharedPtr<Scalar> wf, ErrorCalculator<Scalar>* error_calculator, AdaptivityStoppingCriterion<Scalar>* stopping_criterion_single_step, std::vector<RefinementSelectors::Selector<Scalar>*> selectors, AdaptSolverCriterion* stopping_criterion_global);
      AdaptSolver(SpaceSharedPtr<Scalar> initial_space, WeakFormSharedPtr<Scalar> wf, ErrorCalculator<Scalar>* error_calculator, AdaptivityStoppingCriterion<Scalar>* stopping_criterion_single_step, RefinementSelectors::Selector<Scalar>* selector, AdaptSolverCriterion* stopping_criterion_global);

      void init();

      ~AdaptSolver();

      void solve(AdaptivityType adaptivityType);

      std::vector<MeshFunctionSharedPtr<Scalar> > get_slns();

      MeshFunctionSharedPtr<Scalar> get_sln(int index);

      std::vector<MeshFunctionSharedPtr<Scalar> > get_ref_slns();

      MeshFunctionSharedPtr<Scalar> get_ref_sln(int index);

      void switch_visualization(bool on_off, bool wait_for_keypress);

      void set_exact_solutions(std::vector<MeshFunctionSharedPtr<Scalar> > exact_slns);

      void set_initial_spaces(std::vector<SpaceSharedPtr<Scalar> >);
      void set_wf(WeakFormSharedPtr<Scalar>);
      void set_error_calculator(ErrorCalculator<Scalar>*);
      void set_stopping_criterion_single_step(AdaptivityStoppingCriterion<Scalar>*);
      void set_selectors(std::vector<RefinementSelectors::Selector<Scalar>*>);
      void set_stopping_criterion_global(AdaptSolverCriterion* stopping_criterion_global);

      SolverType* get_solver();
      std::vector<SpaceSharedPtr<Scalar> > get_initial_spaces();
      WeakFormSharedPtr<Scalar> get_wf();
      ErrorCalculator<Scalar>* get_error_calculator();
      AdaptivityStoppingCriterion<Scalar>* get_stopping_criterion_single_step();
      std::vector<RefinementSelectors::Selector<Scalar>*> get_selectors();
      AdaptSolverCriterion* get_stopping_criterion_global();

      virtual void set_verbose_output(bool to_set);

    private:

      virtual bool isOkay() const;
      inline std::string getClassName() const { return "AdaptSolver"; }

      void init_solving(AdaptivityType adaptivityType);

      void deinit_solving();

      void mark_elements_to_reassemble();

      void visualize(std::vector<SpaceSharedPtr<Scalar> >& ref_spaces);

      AdaptSolverCriterion* stopping_criterion_global;

      bool solve_method_running;

      std::vector<SpaceSharedPtr<Scalar> > spaces;

      std::vector<SpaceSharedPtr<Scalar> > ref_spaces;
      std::vector<SpaceSharedPtr<Scalar> > prev_ref_spaces;

      WeakFormSharedPtr<Scalar> wf;

      ErrorCalculator<Scalar>* error_calculator;

      AdaptivityStoppingCriterion<Scalar>* stopping_criterion_single_step;

      std::vector<RefinementSelectors::Selector<Scalar>*> selectors;

      std::vector<MeshFunctionSharedPtr<Scalar> > ref_slns;
      std::vector<MeshFunctionSharedPtr<Scalar> > slns;
      std::vector<MeshFunctionSharedPtr<Scalar> > exact_slns;

      std::vector<Views::ScalarView*> scalar_views;
      std::vector<Views::OrderView*> order_views;
      std::vector<Views::BaseView<Scalar>*> base_views;

      Adapt<Scalar>* adaptivity_internal;

      SolverType* solver;

      unsigned short adaptivity_step;

      std::unordered_set<unsigned int> elements_to_reassemble[H2D_MAX_COMPONENTS];
      std::unordered_set<int> DOFs_to_reassemble[H2D_MAX_COMPONENTS];

      CSCMatrix<Scalar>* prev_mat;
      Vector<Scalar>* prev_rhs, *prev_dirichlet_lift_rhs;

      bool visualization;
      
      bool wait_for_keypress;

      unsigned int total_elements_prev_spaces;

      std::vector<RefinementSelectors::Selector<Scalar>*> hOnlySelectors;
      std::vector<RefinementSelectors::Selector<Scalar>*> pOnlySelectors;

      unsigned char number_of_equations;
    };
  }
}
#endif