problem.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// A generic problem class to keep MH happy
 
// Include guards to prevent multiple inclusion of this header
#ifndef OOMPH_PROBLEM_CLASS_HEADER
#define OOMPH_PROBLEM_CLASS_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// OOMPH-LIB headers
#include "Vector.h"
#include "matrices.h"
#include "generalised_timesteppers.h"
#include "explicit_timesteppers.h"
#include "double_vector_with_halo.h"
#include <complex>
#include <map>
 
 
namespace oomph
{
  // Forward definition for Data class
  class Data;
 
  // Forward definition for Time class
  class Time;
 
  // Forward definition for TimeStepper class
  class TimeStepper;
 
  // Forward definition for Mesh class
  class Mesh;
 
  // Forward definition for RefineableElement class
  class RefineableElement;
 
  // Forward definition for PRefineableElement class
  class PRefineableElement;
 
  // Forward definition for FiniteElement class
  class FiniteElement;
 
  // Forward definition for the GeneralisedElement class
  class GeneralisedElement;
 
  // Forward definition for the Linearsolver class
  class LinearSolver;
 
  // Forward definition for the Eigensolver class
  class EigenSolver;
 
  // Forward definition for the assembly handler
  class AssemblyHandler;
 
  // Forward definition for sum of matrices class
  class SumOfMatrices;
 
 
 
  //=======================================================================
  //=======================================================================
  class Problem : public ExplicitTimeSteppableObject
  {
    // The handler classes need to be friend
    friend class FoldHandler;
    friend class PitchForkHandler;
    friend class HopfHandler;
    template<unsigned NNODE_1D>
    friend class PeriodicOrbitAssemblyHandler;
    friend class BlockFoldLinearSolver;
    friend class BlockPitchForkLinearSolver;
    friend class AugmentedBlockFoldLinearSolver;
    friend class AugmentedBlockPitchForkLinearSolver;
    friend class BlockHopfLinearSolver;
 
 
  private:
    Mesh* Mesh_pt;
 
    Vector<Mesh*> Sub_mesh_pt;
 
    LinearSolver* Linear_solver_pt;
 
    LinearSolver* Mass_matrix_solver_for_explicit_timestepper_pt;
 
    EigenSolver* Eigen_solver_pt;
 
    // Pointer to handler
    AssemblyHandler* Assembly_handler_pt;
 
    LinearSolver* Default_linear_solver_pt;
 
    EigenSolver* Default_eigen_solver_pt;
 
    AssemblyHandler* Default_assembly_handler_pt;
 
    Time* Time_pt;
 
    Vector<TimeStepper*> Time_stepper_pt;
 
    ExplicitTimeStepper* Explicit_time_stepper_pt;
 
    Vector<double>* Saved_dof_pt;
 
    bool Default_set_initial_condition_called;
 
    bool Use_globally_convergent_newton_method;
 
    bool Empty_actions_before_read_unstructured_meshes_has_been_called;
 
    bool Empty_actions_after_read_unstructured_meshes_has_been_called;
 
    bool Store_local_dof_pt_in_elements;
 
    bool Use_predictor_values_as_initial_guess;
 
  protected:
    Vector<Problem*> Copy_of_problem_pt;
 
    std::map<double*, bool> Calculate_dparameter_analytic;
 
    bool Calculate_hessian_products_analytic;
 
  public:
    virtual void debug_hook_fct(const unsigned& i)
    {
      oomph_info << "Called empty hook fct with i=" << i << std::endl;
    }
 
    inline void set_analytic_dparameter(double* const& parameter_pt)
    {
      Calculate_dparameter_analytic[parameter_pt] = true;
    }
 
    inline void unset_analytic_dparameter(double* const& parameter_pt)
    {
      // Find the iterator to the parameter
      std::map<double*, bool>::iterator it =
        Calculate_dparameter_analytic.find(parameter_pt);
      // If the parameter has been found, erase it
      if (it != Calculate_dparameter_analytic.end())
      {
        Calculate_dparameter_analytic.erase(it);
      }
    }
 
    inline bool is_dparameter_calculated_analytically(
      double* const& parameter_pt)
    {
      // If the entry is found then the iterator returned from the find
      // command will NOT be equal to the end and the expression will
      // return true, otherwise it will return false
      return (Calculate_dparameter_analytic.find(parameter_pt) !=
              Calculate_dparameter_analytic.end());
    }
 
    inline void set_analytic_hessian_products()
    {
      Calculate_hessian_products_analytic = true;
    }
 
    void unset_analytic_hessian_products()
    {
      Calculate_hessian_products_analytic = false;
    }
 
    inline bool are_hessian_products_calculated_analytically()
    {
      return Calculate_hessian_products_analytic;
    }
 
    void set_pinned_values_to_zero();
 
 
    static bool Suppress_warning_about_actions_before_read_unstructured_meshes;
 
 
  private:
    double doubly_adaptive_unsteady_newton_solve_helper(
      const double& dt,
      const double& epsilon,
      const unsigned& max_adapt,
      const unsigned& suppress_resolve_after_spatial_adapt,
      const bool& first,
      const bool& shift = true);
 
 
    void refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
                              DocInfo& doc_info,
                              const bool& prune);
 
    void p_refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
                                DocInfo& doc_info,
                                const bool& prune);
 
    void setup_base_mesh_info_after_pruning();
 
    virtual void sparse_assemble_row_or_column_compressed_with_vectors_of_pairs(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    virtual void sparse_assemble_row_or_column_compressed_with_two_vectors(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    virtual void sparse_assemble_row_or_column_compressed_with_maps(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    virtual void sparse_assemble_row_or_column_compressed_with_lists(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    virtual void sparse_assemble_row_or_column_compressed_with_two_arrays(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    Vector<Data*> Global_data_pt;
 
    void calculate_continuation_derivatives_helper(const DoubleVector& z);
 
    void calculate_continuation_derivatives_fd_helper(
      double* const& parameter_pt);
 
    void bifurcation_adapt_helper(unsigned& n_refined,
                                  unsigned& n_unrefined,
                                  const unsigned& bifurcation_type,
                                  const bool& actually_adapt = true);
 
    void bifurcation_adapt_doc_errors(const unsigned& bifurcation_type);
 
    // ALH_TEMP_DEVELOPMENT
  protected:
    LinearAlgebraDistribution* Dof_distribution_pt;
 
  private:
#ifdef OOMPH_HAS_MPI
 
    void get_my_eqns(AssemblyHandler* const& assembly_handler_pt,
                     const unsigned& el_lo,
                     const unsigned& el_hi,
                     Vector<unsigned>& my_eqns);
 
    void parallel_sparse_assemble(
      const LinearAlgebraDistribution* const& dist_pt,
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residuals);
 
    void copy_haloed_eqn_numbers_helper(const bool& do_halos,
                                        const bool& do_external_halos);
 
    void remove_duplicate_data(Mesh* const& mesh_pt,
                               bool& actually_removed_some_data);
 
 
    void remove_null_pointers_from_external_halo_node_storage();
 
    void recompute_load_balanced_assembly();
 
    bool Doc_imbalance_in_parallel_assembly;
 
    bool Use_default_partition_in_load_balance;
 
    Vector<unsigned> First_el_for_assembly;
 
    Vector<unsigned> Last_el_plus_one_for_assembly;
 
    bool Must_recompute_load_balance_for_assembly;
 
    std::map<GeneralisedElement*, unsigned> Base_mesh_element_number_plus_one;
 
 
    Vector<GeneralisedElement*> Base_mesh_element_pt;
 
#endif
 
  protected:
    Vector<double*> Dof_pt;
 
    DoubleVectorWithHaloEntries Element_count_per_dof;
 
    unsigned setup_element_count_per_dof();
 
 
#ifdef OOMPH_HAS_MPI
 
    Vector<double> Elemental_assembly_time;
 
    DoubleVectorHaloScheme* Halo_scheme_pt;
 
    Vector<double*> Halo_dof_pt;
 
    void setup_dof_halo_scheme();
 
#endif
 
    //--------------------- Newton solver parameters
 
    double Relaxation_factor;
 
    double Newton_solver_tolerance;
 
    unsigned Max_newton_iterations;
 
    unsigned Nnewton_iter_taken;
 
    Vector<double> Max_res;
 
    double Max_residuals;
 
    bool Time_adaptive_newton_crash_on_solve_fail;
 
    bool Jacobian_reuse_is_enabled;
 
    bool Jacobian_has_been_computed;
 
    bool Problem_is_nonlinear;
 
    bool Pause_at_end_of_sparse_assembly;
 
    bool Doc_time_in_distribute;
 
    unsigned Sparse_assembly_method;
 
    enum Assembly_method
    {
      Perform_assembly_using_vectors_of_pairs,
      Perform_assembly_using_two_vectors,
      Perform_assembly_using_maps,
      Perform_assembly_using_lists,
      Perform_assembly_using_two_arrays
    };
 
    unsigned Sparse_assemble_with_arrays_initial_allocation;
 
    unsigned Sparse_assemble_with_arrays_allocation_increment;
 
    Vector<Vector<unsigned>> Sparse_assemble_with_arrays_previous_allocation;
 
    double Numerical_zero_for_sparse_assembly;
 
    virtual void sparse_assemble_row_or_column_compressed(
      Vector<int*>& column_or_row_index,
      Vector<int*>& row_or_column_start,
      Vector<double*>& value,
      Vector<unsigned>& nnz,
      Vector<double*>& residual,
      bool compressed_row_flag);
 
    double FD_step_used_in_get_hessian_vector_products;
 
    //---------------------Explicit time-stepping parameters
 
    bool Mass_matrix_reuse_is_enabled;
 
    bool Mass_matrix_has_been_computed;
 
    //---------------------Discontinuous control flags
    bool Discontinuous_element_formulation;
 
 
    //--------------------- Adaptive time-stepping parameters
 
    double Minimum_dt;
 
    double Maximum_dt;
 
    double DTSF_max_increase;
 
    double DTSF_min_decrease;
 
    double Minimum_dt_but_still_proceed;
 
 
    //---------------------  Arc-length continuation paramaters
 
    bool Scale_arc_length;
 
    double Desired_proportion_of_arc_length;
 
    double Theta_squared;
 
    int Sign_of_jacobian;
 
    double Continuation_direction;
 
    double Parameter_derivative;
 
    double Parameter_current;
 
    bool Use_continuation_timestepper;
 
    static ContinuationStorageScheme Continuation_time_stepper;
 
    unsigned Dof_derivative_offset;
 
    unsigned Dof_current_offset;
 
    Vector<double> Dof_derivative;
 
    Vector<double> Dof_current;
 
    double Ds_current;
 
    unsigned Desired_newton_iterations_ds;
 
    double Minimum_ds;
 
    bool Bifurcation_detection;
 
    bool Bisect_to_find_bifurcation;
 
    bool First_jacobian_sign_change;
 
    bool Arc_length_step_taken;
 
    bool Use_finite_differences_for_continuation_derivatives;
 
  public:
    bool distributed() const
    {
#ifdef OOMPH_HAS_MPI
      return Problem_has_been_distributed;
#else
      return false;
#endif
    }
 
#ifdef OOMPH_HAS_MPI
 
 
    enum Distributed_problem_matrix_distribution{Default_matrix_distribution,
                                                 Problem_matrix_distribution,
                                                 Uniform_matrix_distribution};
 
    Distributed_problem_matrix_distribution& distributed_problem_matrix_distribution()
    {
      return Dist_problem_matrix_distribution;
    }
 
    void enable_doc_imbalance_in_parallel_assembly()
    {
      Doc_imbalance_in_parallel_assembly = true;
    }
 
    void disable_doc_imbalance_in_parallel_assembly()
    {
      Doc_imbalance_in_parallel_assembly = false;
    }
 
    Vector<double> elemental_assembly_time()
    {
      return Elemental_assembly_time;
    }
 
    void clear_elemental_assembly_time()
    {
      Must_recompute_load_balance_for_assembly = true;
      Elemental_assembly_time.clear();
    }
 
  private:
    void get_data_to_be_sent_during_load_balancing(
      const Vector<unsigned>& element_domain_on_this_proc,
      Vector<int>& send_n,
      Vector<double>& send_data,
      Vector<int>& send_displacement,
      Vector<unsigned>& old_domain_for_base_element,
      Vector<unsigned>& new_domain_for_base_element,
      unsigned& max_refinement_level_overall);
 
    void get_flat_packed_refinement_pattern_for_load_balancing(
      const Vector<unsigned>& old_domain_for_base_element,
      const Vector<unsigned>& new_domain_for_base_element,
      const unsigned& max_refinement_level_overall,
      std::map<unsigned, Vector<unsigned>>&
        flat_packed_refinement_info_for_root);
 
 
    void send_data_to_be_sent_during_load_balancing(
      Vector<int>& send_n,
      Vector<double>& send_data,
      Vector<int>& send_displacement);
 
    void send_refinement_info_helper(
      Vector<unsigned>& old_domain_for_base_element,
      Vector<unsigned>& new_domain_for_base_element,
      const unsigned& max_refinement_level_overall,
      std::map<unsigned, Vector<unsigned>>& refinement_info_for_root_local,
      Vector<Vector<Vector<unsigned>>>& refinement_info_for_root_elements);
 
    Distributed_problem_matrix_distribution Dist_problem_matrix_distribution;
 
    unsigned Parallel_sparse_assemble_previous_allocation;
 
  protected:
    bool Problem_has_been_distributed;
 
    bool Bypass_increase_in_dof_check_during_pruning;
 
  public:
    void enable_problem_distributed()
    {
      Problem_has_been_distributed = true;
    }
 
    void disable_problem_distributed()
    {
      Problem_has_been_distributed = false;
    }
 
    void set_default_first_and_last_element_for_assembly();
 
    void set_first_and_last_element_for_assembly(
      Vector<unsigned>& first_el_for_assembly,
      Vector<unsigned>& last_el_for_assembly)
    {
      First_el_for_assembly = first_el_for_assembly;
      Last_el_plus_one_for_assembly = last_el_for_assembly;
      unsigned n = last_el_for_assembly.size();
      for (unsigned i = 0; i < n; i++)
      {
        Last_el_plus_one_for_assembly[i] += 1;
      }
    }
 
 
    void get_first_and_last_element_for_assembly(
      Vector<unsigned>& first_el_for_assembly,
      Vector<unsigned>& last_el_for_assembly) const
    {
      first_el_for_assembly = First_el_for_assembly;
      last_el_for_assembly = Last_el_plus_one_for_assembly;
      unsigned n = last_el_for_assembly.size();
      for (unsigned i = 0; i < n; i++)
      {
        last_el_for_assembly[i] -= 1;
      }
    }
 
#endif
 
    virtual void actions_before_adapt() {}
 
    virtual void actions_after_adapt() {}
 
  protected:
    virtual void actions_before_newton_solve() {}
 
    virtual void actions_after_newton_solve() {}
 
    virtual void actions_before_newton_convergence_check() {}
 
    virtual void actions_before_newton_step() {}
 
    virtual void actions_after_newton_step() {}
 
    virtual void actions_before_implicit_timestep() {}
 
    virtual void actions_after_implicit_timestep() {}
 
    virtual void actions_after_implicit_timestep_and_error_estimation() {}
 
    virtual void actions_before_explicit_timestep() {}
 
    virtual void actions_after_explicit_timestep() {}
 
    virtual void actions_before_read_unstructured_meshes()
    {
      Empty_actions_before_read_unstructured_meshes_has_been_called = true;
    }
 
    virtual void actions_after_read_unstructured_meshes()
    {
      Empty_actions_after_read_unstructured_meshes_has_been_called = true;
    }
 
#ifdef OOMPH_HAS_MPI
    virtual void actions_before_distribute() {}
 
    virtual void actions_after_distribute() {}
 
#endif
 
    virtual void actions_after_change_in_global_parameter(
      double* const& parameter_pt)
    {
      // Default action should cover all possibilities
      actions_before_newton_solve();
      actions_before_newton_convergence_check();
      actions_after_newton_solve();
    }
 
    virtual void actions_after_change_in_bifurcation_parameter()
    {
      // Default action should cover all possibilities
      actions_before_newton_solve();
      actions_before_newton_convergence_check();
      actions_after_newton_solve();
    }
 
    virtual void actions_after_parameter_increase(double* const& parameter_pt)
    {
    }
 
    inline double& dof_derivative(const unsigned& i)
    {
      if (!Use_continuation_timestepper)
      {
        return Dof_derivative[i];
      }
      else
      {
        return (*(Dof_pt[i] + Dof_derivative_offset));
      }
    }
 
    inline double& dof_current(const unsigned& i)
    {
      if (!Use_continuation_timestepper)
      {
        return Dof_current[i];
      }
      else
      {
        return (*(Dof_pt[i] + Dof_current_offset));
      }
    }
 
    virtual void set_initial_condition()
    {
      std::ostringstream warn_message;
      warn_message
        << "Warning: We're using the default (empty) set_initial_condition().\n"
        << "If the initial conditions isn't re-assigned after a mesh adaption "
           "\n"
        << "the initial conditions will represent the interpolation of the \n"
        << "initial conditions that were assigned on the original coarse "
           "mesh.\n";
      OomphLibWarning(warn_message.str(),
                      "Problem::set_initial_condition()",
                      OOMPH_EXCEPTION_LOCATION);
 
      // Indicate that this function has been called. This flag is set so
      // that the unsteady_newton_solve routine can be instructed whether
      // or not to shift the history values. If set_initial_condition() has
      // been overloaded than this (default) function won't be called, and
      // so this flag will remain false (its default value). If
      // set_initial_condition() has not been overloaded then this function
      // will be called and the flag set to true.
      Default_set_initial_condition_called = true;
    }
 
    virtual double global_temporal_error_norm()
    {
      std::string error_message =
        "The global_temporal_error_norm function will be problem-specific:\n";
      error_message += "Please write your own in your Problem class";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      return 0.0;
    }
 
    OomphCommunicator* Communicator_pt;
 
  public:
    OomphCommunicator* communicator_pt()
    {
      return Communicator_pt;
    }
 
    const OomphCommunicator* communicator_pt() const
    {
      return Communicator_pt;
    }
 
    typedef void (*SpatialErrorEstimatorFctPt)(Mesh*& mesh_pt,
                                               Vector<double>& elemental_error);
 
    typedef void (*SpatialErrorEstimatorWithDocFctPt)(
      Mesh*& mesh_pt, Vector<double>& elemental_error, DocInfo& doc_info);
 
    Problem();
 
    Problem(const Problem& dummy) = delete;
 
    void operator=(const Problem&) = delete;
 
    virtual ~Problem();
 
    Mesh*& mesh_pt()
    {
      return Mesh_pt;
    }
 
    Mesh* const& mesh_pt() const
    {
      return Mesh_pt;
    }
 
    Mesh*& mesh_pt(const unsigned& imesh)
    {
      // If there are no submeshes then the zero-th submesh is the
      // mesh itself
      if ((Sub_mesh_pt.size() == 0) && (imesh == 0))
      {
        return Mesh_pt;
      }
      else
      {
        return Sub_mesh_pt[imesh];
      }
    }
 
    Mesh* const& mesh_pt(const unsigned& imesh) const
    {
      // If there are no submeshes then the zero-th submesh is the
      // mesh itself
      if ((Sub_mesh_pt.size() == 0) && (imesh == 0))
      {
        return Mesh_pt;
      }
      else
      {
        return Sub_mesh_pt[imesh];
      }
    }
 
    unsigned nsub_mesh() const
    {
      return Sub_mesh_pt.size();
    }
 
    unsigned add_sub_mesh(Mesh* const& mesh_pt)
    {
      Sub_mesh_pt.push_back(mesh_pt);
      return Sub_mesh_pt.size() - 1;
    }
 
 
    void flush_sub_meshes()
    {
      Sub_mesh_pt.clear();
    }
 
    void build_global_mesh();
 
    void rebuild_global_mesh();
 
#ifdef OOMPH_HAS_MPI
 
    virtual void build_mesh()
    {
      std::string error_message = "You are likely to have reached this error "
                                  "from Problem::load_balance()\n";
      error_message +=
        "The build_mesh() function is problem-specific and must be supplied:\n";
      error_message +=
        "Please write your own in your Problem class, ensuring that\n";
      error_message += "any (sub)meshes are built before calling "
                       "Problem::build_global_mesh().\n";
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void load_balance()
    {
      // Dummy DocInfo
      DocInfo doc_info;
      doc_info.disable_doc();
 
      // Don't report stats
      bool report_stats = false;
 
      // Dummy imposed partioning vector
      Vector<unsigned> input_target_domain_for_local_non_halo_element;
 
      // Do it!
      load_balance(
        doc_info, report_stats, input_target_domain_for_local_non_halo_element);
    }
 
    void load_balance(const bool& report_stats)
    {
      // Dummy DocInfo
      DocInfo doc_info;
      doc_info.disable_doc();
 
      // Dummy imposed partioning vector
      Vector<unsigned> input_target_domain_for_local_non_halo_element;
 
      // Do it
      load_balance(
        doc_info, report_stats, input_target_domain_for_local_non_halo_element);
    }
 
 
    void load_balance(DocInfo& doc_info, const bool& report_stats)
    {
      // Dummy imposed partioning vector
      Vector<unsigned> input_target_domain_for_local_non_halo_element;
 
      // Do it
      load_balance(
        doc_info, report_stats, input_target_domain_for_local_non_halo_element);
    }
 
    void load_balance(
      DocInfo& doc_info,
      const bool& report_stats,
      const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
 
    void set_default_partition_in_load_balance()
    {
      Use_default_partition_in_load_balance = true;
    }
 
    void unset_default_partition_in_load_balance()
    {
      Use_default_partition_in_load_balance = false;
    }
 
    void refine_distributed_base_mesh(
      Vector<Vector<Vector<unsigned>>>& to_be_refined_on_each_root,
      const unsigned& max_level_overall);
 
#endif
 
    LinearSolver*& linear_solver_pt()
    {
      return Linear_solver_pt;
    }
 
    LinearSolver* const& linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
 
    LinearSolver*& mass_matrix_solver_for_explicit_timestepper_pt()
    {
      return Mass_matrix_solver_for_explicit_timestepper_pt;
    }
 
    LinearSolver* mass_matrix_solver_for_explicit_timestepper_pt() const
    {
      return Mass_matrix_solver_for_explicit_timestepper_pt;
    }
 
    EigenSolver*& eigen_solver_pt()
    {
      return Eigen_solver_pt;
    }
 
    EigenSolver* const& eigen_solver_pt() const
    {
      return Eigen_solver_pt;
    }
 
    Time*& time_pt()
    {
      return Time_pt;
    }
 
    Time* time_pt() const
    {
      return Time_pt;
    }
 
    double& time();
 
    double time() const;
 
 
    TimeStepper*& time_stepper_pt()
    {
      if (Time_stepper_pt.size() == 0)
      {
        throw OomphLibError("No timestepper allocated yet\n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return Time_stepper_pt[0];
    }
 
    const TimeStepper* time_stepper_pt() const
    {
      if (Time_stepper_pt.size() == 0)
      {
        throw OomphLibError("No timestepper allocated yet\n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return Time_stepper_pt[0];
    }
 
    TimeStepper*& time_stepper_pt(const unsigned& i)
    {
      return Time_stepper_pt[i];
    }
 
    ExplicitTimeStepper*& explicit_time_stepper_pt()
    {
      return Explicit_time_stepper_pt;
    }
 
    unsigned long set_timestepper_for_all_data(
      TimeStepper* const& time_stepper_pt,
      const bool& preserve_existing_data = false);
 
    virtual void shift_time_values();
 
    AssemblyHandler*& assembly_handler_pt()
    {
      return Assembly_handler_pt;
    }
 
    AssemblyHandler* const& assembly_handler_pt() const
    {
      return Assembly_handler_pt;
    }
 
    double& minimum_dt()
    {
      return Minimum_dt;
    }
 
    double& maximum_dt()
    {
      return Maximum_dt;
    }
 
    unsigned& max_newton_iterations()
    {
      return Max_newton_iterations;
    }
 
    void problem_is_nonlinear(const bool& prob_lin)
    {
      Problem_is_nonlinear = prob_lin;
    }
 
 
    double& max_residuals()
    {
      return Max_residuals;
    }
 
    bool& time_adaptive_newton_crash_on_solve_fail()
    {
      return Time_adaptive_newton_crash_on_solve_fail;
    }
 
    double& newton_solver_tolerance()
    {
      return Newton_solver_tolerance;
    }
 
    void add_time_stepper_pt(TimeStepper* const& time_stepper_pt);
 
    void set_explicit_time_stepper_pt(
      ExplicitTimeStepper* const& explicit_time_stepper_pt);
 
 
    void initialise_dt(const double& dt);
 
    void initialise_dt(const Vector<double>& dt);
 
    Data*& global_data_pt(const unsigned& i)
    {
      return Global_data_pt[i];
    }
 
    void add_global_data(Data* const& global_data_pt)
    {
      Global_data_pt.push_back(global_data_pt);
    }
 
 
    void flush_global_data()
    {
      Global_data_pt.resize(0);
    }
 
    LinearAlgebraDistribution* const& dof_distribution_pt() const
    {
      return Dof_distribution_pt;
    }
 
    unsigned long ndof() const
    {
      return Dof_distribution_pt->nrow();
    }
 
    unsigned ntime_stepper() const
    {
      return Time_stepper_pt.size();
    }
 
    unsigned nglobal_data() const
    {
      return Global_data_pt.size();
    }
 
    unsigned self_test();
 
    void enable_store_local_dof_pt_in_elements()
    {
      Store_local_dof_pt_in_elements = true;
    }
 
    void disable_store_local_dof_pt_in_elements()
    {
      Store_local_dof_pt_in_elements = false;
    }
 
    unsigned long assign_eqn_numbers(
      const bool& assign_local_eqn_numbers = true);
 
    void describe_dofs(std::ostream& out = *(oomph_info.stream_pt())) const;
 
    void enable_discontinuous_formulation()
    {
      Discontinuous_element_formulation = true;
    }
 
    void disable_discontinuous_formulation()
    {
      Discontinuous_element_formulation = false;
    }
 
    void get_dofs(DoubleVector& dofs) const;
 
    void get_dofs(const unsigned& t, DoubleVector& dofs) const;
 
    void set_dofs(const DoubleVector& dofs);
 
    void set_dofs(const unsigned& t, DoubleVector& dofs);
 
    void set_dofs(const unsigned& t, Vector<double*>& dof_pt);
 
    void add_to_dofs(const double& lambda, const DoubleVector& increment_dofs);
 
    inline double* global_dof_pt(const unsigned& i)
    {
#ifdef OOMPH_HAS_MPI
      // If the problem is distributed I have to do something different
      if (Problem_has_been_distributed)
      {
#ifdef PARANOID
        if (Halo_scheme_pt == 0)
        {
          std::ostringstream error_stream;
          error_stream
            << "Direct access to global dofs in a distributed problem is only\n"
            << "possible if the function setup_dof_halo_scheme() has been "
               "called.\n"
            << "You can access local dofs via Problem::dof(i).\n";
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Work out the local coordinate of the dof
        // based on the original distribution stored in the Halo_scheme
        const unsigned first_row_local =
          Halo_scheme_pt->distribution_pt()->first_row();
        const unsigned n_row_local =
          Halo_scheme_pt->distribution_pt()->nrow_local();
 
        // If we are in range then just call the local value
        if ((i >= first_row_local) && (i < first_row_local + n_row_local))
        {
          return this->Dof_pt[i - first_row_local];
        }
        // Otherwise the entry is not stored in the local processor
        // and we must have haloed it
        else
        {
          return Halo_dof_pt[Halo_scheme_pt->local_index(i)];
        }
      }
      // Otherwise just return the dof
      else
#endif
      {
        return this->Dof_pt[i];
      }
    }
 
    double& dof(const unsigned& i)
    {
      return *(Dof_pt[i]);
    }
 
    double dof(const unsigned& i) const
    {
      return *(Dof_pt[i]);
    }
 
    double*& dof_pt(const unsigned& i)
    {
      return Dof_pt[i];
    }
 
    double* dof_pt(const unsigned& i) const
    {
      return Dof_pt[i];
    }
 
    virtual void get_inverse_mass_matrix_times_residuals(DoubleVector& Mres);
 
    virtual void get_dvaluesdt(DoubleVector& f);
 
    virtual void get_residuals(DoubleVector& residuals);
 
    virtual void get_jacobian(DoubleVector& residuals,
                              DenseDoubleMatrix& jacobian);
 
    virtual void get_jacobian(DoubleVector& residuals,
                              CRDoubleMatrix& jacobian);
 
    virtual void get_jacobian(DoubleVector& residuals,
                              CCDoubleMatrix& jacobian);
 
    virtual void get_jacobian(DoubleVector& residuals, SumOfMatrices& jacobian)
    {
      std::ostringstream error_stream;
      error_stream
        << "You must overload this function in your problem class to specify\n"
        << "which matrices should be summed to give the final Jacobian.";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void get_fd_jacobian(DoubleVector& residuals,
                         DenseMatrix<double>& jacobian);
 
    void get_derivative_wrt_global_parameter(double* const& parameter_pt,
                                             DoubleVector& result);
 
    void get_hessian_vector_products(
      DoubleVectorWithHaloEntries const& Y,
      Vector<DoubleVectorWithHaloEntries> const& C,
      Vector<DoubleVectorWithHaloEntries>& product);
 
 
    // void get_derivative_wrt_global_parameter(double* const &parameter_pt,
    //                                         GeneralisedElement* const
    //                                         &elem_pt, Vector<double>
    //                                         &result);
 
    void solve_eigenproblem(const unsigned& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector,
                            const bool& steady = true);
 
    void solve_eigenproblem(const unsigned& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            const bool& steady = true)
    {
      // Create temporary storage for the eigenvectors (potentially wasteful)
      Vector<DoubleVector> eigenvector;
      solve_eigenproblem(n_eval, eigenvalue, eigenvector, steady);
    }
 
    virtual void get_eigenproblem_matrices(CRDoubleMatrix& mass_matrix,
                                           CRDoubleMatrix& main_matrix,
                                           const double& shift = 0.0);
 
    void assign_eigenvector_to_dofs(DoubleVector& eigenvector);
 
    void add_eigenvector_to_dofs(const double& epsilon,
                                 const DoubleVector& eigenvector);
 
 
    void store_current_dof_values();
 
    void restore_dof_values();
 
    void enable_jacobian_reuse()
    {
      Jacobian_reuse_is_enabled = true;
      Jacobian_has_been_computed = false;
    }
 
    void disable_jacobian_reuse()
    {
      Jacobian_reuse_is_enabled = false;
      Jacobian_has_been_computed = false;
    }
 
    bool jacobian_reuse_is_enabled()
    {
      return Jacobian_reuse_is_enabled;
    }
 
    bool& use_predictor_values_as_initial_guess()
    {
      return Use_predictor_values_as_initial_guess;
    }
 
    void newton_solve();
 
    void enable_globally_convergent_newton_method()
    {
      Use_globally_convergent_newton_method = true;
    }
 
    void disable_globally_convergent_newton_method()
    {
      Use_globally_convergent_newton_method = false;
    }
 
  private:
    void globally_convergent_line_search(
      const Vector<double>& x_old,
      const double& half_residual_squared_old,
      DoubleVector& gradient,
      DoubleVector& newton_dir,
      double& half_residual_squared,
      const double& stpmax);
 
  public:
    void newton_solve(unsigned const& max_adapt);
 
    void steady_newton_solve(unsigned const& max_adapt = 0);
 
    void copy(Problem* orig_problem_pt);
 
    virtual Problem* make_copy();
 
    virtual void read(std::ifstream& restart_file, bool& unsteady_restart);
 
    virtual void read(std::ifstream& restart_file)
    {
      bool unsteady_restart;
      read(restart_file, unsteady_restart);
    }
 
    virtual void dump(std::ofstream& dump_file) const;
 
    void dump(const std::string& dump_file_name) const
    {
      std::ofstream dump_stream(dump_file_name.c_str());
#ifdef PARANOID
      if (!dump_stream.is_open())
      {
        std::string err = "Couldn't open file " + dump_file_name;
        throw OomphLibError(
          err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
      }
#endif
      dump(dump_stream);
    }
 
#ifdef OOMPH_HAS_MPI
 
    void get_all_halo_data(std::map<unsigned, double*>& map_of_halo_data);
 
    long synchronise_eqn_numbers(const bool& assign_local_eqn_numbers = true);
 
    void synchronise_dofs(const bool& do_halos, const bool& do_external_halos);
 
    void synchronise_all_dofs();
 
    void check_halo_schemes(DocInfo& doc_info);
 
    void check_halo_schemes()
    {
      DocInfo tmp_doc_info;
      tmp_doc_info.disable_doc();
      check_halo_schemes(tmp_doc_info);
    }
 
    Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
                                DocInfo& doc_info,
                                const bool& report_stats = false);
 
    Vector<unsigned> distribute(DocInfo& doc_info,
                                const bool& report_stats = false);
 
    Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
                                const bool& report_stats = false);
 
    Vector<unsigned> distribute(const bool& report_stats = false);
 
    virtual void partition_global_mesh(Mesh*& global_mesh_pt,
                                       DocInfo& doc_info,
                                       Vector<unsigned>& element_domain,
                                       const bool& report_stats = false);
 
    void prune_halo_elements_and_nodes(DocInfo& doc_info,
                                       const bool& report_stats);
 
    void prune_halo_elements_and_nodes(const bool& report_stats = false)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      prune_halo_elements_and_nodes(doc_info, report_stats);
    }
 
    double Max_permitted_error_for_halo_check;
 
    bool problem_has_been_distributed()
    {
      return Problem_has_been_distributed;
    }
 
#endif
 
    void delete_all_external_storage();
 
    bool Shut_up_in_newton_solve;
 
 
  protected:
    bool Always_take_one_newton_step;
 
    double Timestep_reduction_factor_after_nonconvergence;
 
 
    bool Keep_temporal_error_below_tolerance;
 
    unsigned newton_solve_continuation(double* const& parameter_pt);
 
    unsigned newton_solve_continuation(double* const& parameter_pt,
                                       DoubleVector& z);
 
    void calculate_continuation_derivatives(double* const& parameter_pt);
 
    void calculate_continuation_derivatives(const DoubleVector& z);
 
    void calculate_continuation_derivatives_fd(double* const& parameter_pt);
 
    bool does_pointer_correspond_to_problem_data(double* const& parameter_pt);
 
  public:
    virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking();
 
    double* bifurcation_parameter_pt() const;
 
 
    void get_bifurcation_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
 
    void activate_fold_tracking(double* const& parameter_pt,
                                const bool& block_solve = true);
 
    void activate_bifurcation_tracking(double* const& parameter_pt,
                                       const DoubleVector& eigenvector,
                                       const bool& block_solve = true);
 
    void activate_bifurcation_tracking(double* const& parameter_pt,
                                       const DoubleVector& eigenvector,
                                       const DoubleVector& normalisation,
                                       const bool& block_solve = true);
 
 
    void activate_pitchfork_tracking(double* const& parameter_pt,
                                     const DoubleVector& symmetry_vector,
                                     const bool& block_solve = true);
 
    void activate_hopf_tracking(double* const& parameter_pt,
                                const bool& block_solve = true);
 
    void activate_hopf_tracking(double* const& parameter_pt,
                                const double& omega,
                                const DoubleVector& null_real,
                                const DoubleVector& null_imag,
                                const bool& block_solve = true);
 
    void deactivate_bifurcation_tracking()
    {
      reset_assembly_handler_to_default();
    }
 
    void reset_assembly_handler_to_default();
 
  private:
    double arc_length_step_solve_helper(double* const& parameter_pt,
                                        const double& ds,
                                        const unsigned& max_adapt);
 
 
    // ALH_DEVELOP
  protected:
    void set_consistent_pinned_values_for_continuation();
 
  public:
    double arc_length_step_solve(double* const& parameter_pt,
                                 const double& ds,
                                 const unsigned& max_adapt = 0);
 
    double arc_length_step_solve(Data* const& data_pt,
                                 const unsigned& data_index,
                                 const double& ds,
                                 const unsigned& max_adapt = 0);
 
 
    void reset_arc_length_parameters()
    {
      Theta_squared = 1.0;
      Sign_of_jacobian = 0;
      Continuation_direction = 1.0;
      Parameter_derivative = 1.0;
      First_jacobian_sign_change = false;
      Arc_length_step_taken = false;
      Dof_derivative.resize(0);
    }
 
    int& sign_of_jacobian()
    {
      return Sign_of_jacobian;
    }
 
    void explicit_timestep(const double& dt, const bool& shift_values = true);
 
    void unsteady_newton_solve(const double& dt);
 
    void unsteady_newton_solve(const double& dt, const bool& shift_values);
 
    void unsteady_newton_solve(const double& dt,
                               const unsigned& max_adapt,
                               const bool& first,
                               const bool& shift = true);
 
 
    double doubly_adaptive_unsteady_newton_solve(const double& dt,
                                                 const double& epsilon,
                                                 const unsigned& max_adapt,
                                                 const bool& first,
                                                 const bool& shift = true)
    {
      // Call helper function with default setting (do re-solve after
      // spatial adaptation)
      unsigned suppress_resolve_after_spatial_adapt_flag = 0;
      return doubly_adaptive_unsteady_newton_solve_helper(
        dt,
        epsilon,
        max_adapt,
        suppress_resolve_after_spatial_adapt_flag,
        first,
        shift);
    }
 
 
    double doubly_adaptive_unsteady_newton_solve(
      const double& dt,
      const double& epsilon,
      const unsigned& max_adapt,
      const unsigned& suppress_resolve_after_spatial_adapt_flag,
      const bool& first,
      const bool& shift = true)
    {
      // Call helper function
      return doubly_adaptive_unsteady_newton_solve_helper(
        dt,
        epsilon,
        max_adapt,
        suppress_resolve_after_spatial_adapt_flag,
        first,
        shift);
    }
 
 
    double adaptive_unsteady_newton_solve(const double& dt_desired,
                                          const double& epsilon);
 
    double adaptive_unsteady_newton_solve(const double& dt_desired,
                                          const double& epsilon,
                                          const bool& shift_values);
 
    void assign_initial_values_impulsive();
 
    void assign_initial_values_impulsive(const double& dt);
 
    void calculate_predictions();
 
    void enable_mass_matrix_reuse();
 
    void disable_mass_matrix_reuse();
 
    bool mass_matrix_reuse_is_enabled()
    {
      return Mass_matrix_reuse_is_enabled;
    }
 
    void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      bool prune = false;
      refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
                          DocInfo& doc_info)
    {
      bool prune = false;
      refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      bool prune = true;
      refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh,
                                    DocInfo& doc_info)
    {
      bool prune = true;
      refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void refine_uniformly(DocInfo& doc_info)
    {
      // Number of (sub)meshes
      unsigned nmesh = std::max(unsigned(1), nsub_mesh());
 
      // Refine each mesh once
      Vector<unsigned> nrefine_for_mesh(nmesh, 1);
      refine_uniformly(nrefine_for_mesh);
    }
 
 
    void refine_uniformly_and_prune(DocInfo& doc_info)
    {
      // Number of (sub)meshes
      unsigned nmesh = std::max(unsigned(1), nsub_mesh());
 
      // Refine each mesh once
      Vector<unsigned> nrefine_for_mesh(nmesh, 1);
      refine_uniformly_and_prune(nrefine_for_mesh);
    }
 
 
    void refine_uniformly()
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      refine_uniformly(doc_info);
    }
 
    void refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
 
    void refine_uniformly(const unsigned& i_mesh)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      refine_uniformly(i_mesh, doc_info);
    }
 
    void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      bool prune = false;
      p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
                            DocInfo& doc_info)
    {
      bool prune = false;
      p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void p_refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh)
    {
      // Not tested:
      throw OomphLibWarning("This functionality has not yet been tested.",
                            "Problem::p_refine_uniformly_and_prune()",
                            OOMPH_EXCEPTION_LOCATION);
      DocInfo doc_info;
      doc_info.disable_doc();
      bool prune = true;
      p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void p_refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh,
                                      DocInfo& doc_info)
    {
      // Not tested:
      throw OomphLibWarning("This functionality has not yet been tested.",
                            "Problem::p_refine_uniformly_and_prune()",
                            OOMPH_EXCEPTION_LOCATION);
      bool prune = true;
      p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
    }
 
    void p_refine_uniformly(DocInfo& doc_info)
    {
      // Number of (sub)meshes
      unsigned nmesh = std::max(unsigned(1), nsub_mesh());
 
      // Refine each mesh once
      Vector<unsigned> nrefine_for_mesh(nmesh, 1);
      p_refine_uniformly(nrefine_for_mesh);
    }
 
 
    void p_refine_uniformly_and_prune(DocInfo& doc_info)
    {
      // Not tested:
      throw OomphLibWarning("This functionality has not yet been tested.",
                            "Problem::p_refine_uniformly_and_prune()",
                            OOMPH_EXCEPTION_LOCATION);
      // Number of (sub)meshes
      unsigned nmesh = std::max(unsigned(1), nsub_mesh());
 
      // Refine each mesh once
      Vector<unsigned> nrefine_for_mesh(nmesh, 1);
      p_refine_uniformly_and_prune(nrefine_for_mesh);
    }
 
 
    void p_refine_uniformly()
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      p_refine_uniformly(doc_info);
    }
 
    void p_refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
 
    void p_refine_uniformly(const unsigned& i_mesh)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      p_refine_uniformly(i_mesh, doc_info);
    }
 
    void refine_selected_elements(
      const Vector<unsigned>& elements_to_be_refined);
 
 
    void refine_selected_elements(
      const Vector<RefineableElement*>& elements_to_be_refined_pt);
 
    void refine_selected_elements(
      const unsigned& i_mesh, const Vector<unsigned>& elements_to_be_refined);
 
    void refine_selected_elements(
      const unsigned& i_mesh,
      const Vector<RefineableElement*>& elements_to_be_refined_pt);
 
    void refine_selected_elements(
      const Vector<Vector<unsigned>>& elements_to_be_refined);
 
    void refine_selected_elements(
      const Vector<Vector<RefineableElement*>>& elements_to_be_refined_pt);
 
    void p_refine_selected_elements(
      const Vector<unsigned>& elements_to_be_refined);
 
 
    void p_refine_selected_elements(
      const Vector<PRefineableElement*>& elements_to_be_refined_pt);
 
    void p_refine_selected_elements(
      const unsigned& i_mesh, const Vector<unsigned>& elements_to_be_refined);
 
    void p_refine_selected_elements(
      const unsigned& i_mesh,
      const Vector<PRefineableElement*>& elements_to_be_refined_pt);
 
    void p_refine_selected_elements(
      const Vector<Vector<unsigned>>& elements_to_be_refined);
 
    void p_refine_selected_elements(
      const Vector<Vector<PRefineableElement*>>& elements_to_be_refined_pt);
 
    unsigned unrefine_uniformly();
 
    unsigned unrefine_uniformly(const unsigned& i_mesh);
 
    void p_unrefine_uniformly(DocInfo& doc_info);
 
    void p_unrefine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
 
    void adapt(unsigned& n_refined, unsigned& n_unrefined);
 
    void adapt()
    {
      unsigned n_refined, n_unrefined;
      adapt(n_refined, n_unrefined);
    }
 
    void p_adapt(unsigned& n_refined, unsigned& n_unrefined);
 
    void p_adapt()
    {
      unsigned n_refined, n_unrefined;
      p_adapt(n_refined, n_unrefined);
    }
 
 
    void adapt_based_on_error_estimates(
      unsigned& n_refined,
      unsigned& n_unrefined,
      Vector<Vector<double>>& elemental_error);
 
    void adapt_based_on_error_estimates(Vector<Vector<double>>& elemental_error)
    {
      unsigned n_refined, n_unrefined;
      adapt_based_on_error_estimates(n_refined, n_unrefined, elemental_error);
    }
 
 
    void get_all_error_estimates(Vector<Vector<double>>& elemental_error);
 
    void doc_errors(DocInfo& doc_info);
 
    void doc_errors()
    {
      DocInfo tmp_doc_info;
      tmp_doc_info.disable_doc();
      doc_errors(tmp_doc_info);
    }
 
    void enable_info_in_newton_solve()
    {
      Shut_up_in_newton_solve = false;
    }
 
    void disable_info_in_newton_solve()
    {
      Shut_up_in_newton_solve = true;
    }
  };
 
 
  //=======================================================================
  //=======================================================================
  class NewtonSolverError
  {
  public:
    bool linear_solver_error;
 
    unsigned iterations;
 
    double maxres;
 
    NewtonSolverError() : linear_solver_error(false), iterations(0), maxres(0.0)
    {
    }
 
    NewtonSolverError(const bool& Passed_linear_failure)
      : linear_solver_error(Passed_linear_failure), iterations(0), maxres(0.0)
    {
    }
 
    NewtonSolverError(unsigned Passed_iterations, double Passed_maxres)
      : linear_solver_error(false),
        iterations(Passed_iterations),
        maxres(Passed_maxres)
    {
    }
 
    /*  /// Broken copy constructor */
    /*  NewtonSolverError(const NewtonSolverError& dummy) = delete;  */
 
    /*  /// Broken assignment operator */
    /*  void operator=(const NewtonSolverError&) = delete; */
  };
 
 
} // namespace oomph
 
 
#endif