oomph_utilities.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 header containing useful utility classes, functions and constants
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_UTILITIES_HEADER
#define OOMPH_UTILITIES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
#ifdef OOMPH_HAS_MPI
// mpi headers
#include "mpi.h"
#endif
 
// Standard libray headers
#include <string>
#include <sstream>
#include <fstream>
#include <cmath>
#include <map>
#include <ctime>
#include <complex>
 
// oomph-lib headers
#include "Vector.h"
#include "oomph_definitions.h"
#include "communicator.h"
 
/* // Header for cxx-prettyprint: A C++ Container Pretty-Printer. If the */
/* // standard in use is older than C++11 then use a special older version of */
/* // the header. Header is from here: */
/* // http://louisdx.github.io/cxx-prettyprint/ */
/* #if __cplusplus <= 199711L */
/*   #include "prettyprint98.h" */
/* #else */
/*   #include "prettyprint.h" */
/* #endif */
 
// Header for name demangling
#include <cxxabi.h>
 
namespace oomph
{
  //======start_of_ANSIEscapeCode_namespace=================
  //====================================================================
  namespace ANSIEscapeCode
  {
    extern void set_text_effect(std::string text_effect);
 
    extern std::string Text_effect;
 
    extern std::string Black;
    extern std::string Red;
    extern std::string Green;
    extern std::string Yellow;
    extern std::string Blue;
    extern std::string Magenta;
    extern std::string Cyan;
    extern std::string Reset;
  } // namespace ANSIEscapeCode
 
  //=====================================================================
  //=====================================================================
  namespace DebugHelpers
  {
    extern std::string debug_string(const std::string& filename,
                                    const int& line,
                                    const bool& make_new_line = false);
  } // namespace DebugHelpers
 
  // Forward declaration needed for second invariant helper
  template<class TYPE>
  class DenseMatrix;
 
  //=====================================================================
  //=====================================================================
  namespace SecondInvariantHelper
  {
    extern double second_invariant(const DenseMatrix<double>& tensor);
 
  } // namespace SecondInvariantHelper
 
 
  //==============================================
  //==============================================
  namespace BrokenCopy
  {
    extern void broken_assign(const std::string& class_name);
 
    extern void broken_copy(const std::string& class_name);
 
  } // namespace BrokenCopy
 
  //========================================
  //=======================================
  namespace MathematicalConstants
  {
    const double Pi = 3.1415926535897932384626433832795028841971693993751;
 
    const std::complex<double> I(0.0, 1.0);
  } // namespace MathematicalConstants
 
 
  //================================================================
  //================================================================
  template<class T>
  class AbsCmp
  {
  public:
    bool operator()(const T& x, const T& y) const
    {
      return std::abs(x) < std::abs(y);
    }
  };
 
 
  // =================================================================
  // =================================================================
  namespace StringConversion
  {
    template<class T>
    std::string to_string(T object, unsigned float_precision = 8)
    {
      std::stringstream ss;
      ss.precision(float_precision);
      ss << object;
      return ss.str();
    }
 
    std::string to_lower(const std::string& input);
 
    std::string to_upper(const std::string& input);
 
    void split_string(const std::string& s,
                      char delim,
                      Vector<std::string>& elems);
 
    Vector<std::string> split_string(const std::string& s, char delim);
 
  } // namespace StringConversion
 
 
  namespace TypeNames
  {
    template<class T>
    std::string get_type_name(T& obj)
    {
      using namespace StringConversion;
 
      int status;
      std::string typestr =
        to_string(abi::__cxa_demangle(typeid(obj).name(), 0, 0, &status));
 
      // Throw if there was an error
      if (status != 0)
      {
        std::string err = "Error code " + to_string(status) +
                          " in demangler, see documentation of "
                          "abi::__cxa_demangle for the meaning";
        throw OomphLibError(
          err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      return typestr;
    }
 
    template<class T>
    std::string get_type_name(T* obj)
    {
      return get_type_name(*obj);
    }
 
  } // namespace TypeNames
 
 
 
  //====================================================================
  //====================================================================
  namespace CumulativeTimings
  {
    extern void start(const unsigned& i);
 
    extern void halt(const unsigned& i);
 
    extern void reset(const unsigned& i);
 
    extern void reset();
 
    extern double cumulative_time(const unsigned& i);
 
    extern void set_ntimers(const unsigned& ntimers);
 
    extern Vector<clock_t> Timing;
 
    extern Vector<clock_t> Start_time;
 
  } // namespace CumulativeTimings
 
 
  // ============================================================
  // Automatically checked casting functions (from boost)
  // ============================================================
 
  template<class Target, class Source>
  inline Target checked_dynamic_cast(Source* x)
  {
    Target tmp = dynamic_cast<Target>(x);
#ifdef PARANOID
    if (tmp == 0)
    {
      throw OomphLibError(
        "Bad cast.", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
    return tmp;
  }
 
  template<class Target, class Source>
  inline Target checked_static_cast(Source* x)
  {
#ifdef PARANOID
    // Check that the cast will work as expected.
    if (dynamic_cast<Target>(x) != x)
    {
      throw OomphLibError(
        "Bad cast.", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
    return static_cast<Target>(x);
  }
 
 
  //====================================================================
  //====================================================================
  class Timer
  {
  public:
    Timer(const unsigned& n_timer)
    {
      set_ntimers(n_timer);
    }
 
    void start(const unsigned& i)
    {
      Start_time[i] = clock();
    }
 
    void halt(const unsigned& i)
    {
      Timing[i] += clock() - Start_time[i];
    }
 
    double cumulative_time(const unsigned& i)
    {
      return double(Timing[i]) / CLOCKS_PER_SEC;
    }
 
    void reset(const unsigned& i)
    {
      Timing[i] = clock_t(0.0);
    }
 
    void reset()
    {
      unsigned n = Timing.size();
      for (unsigned i = 0; i < n; i++)
      {
        Timing[i] = clock_t(0.0);
      }
    }
 
    void set_ntimers(const unsigned& ntimers)
    {
      Timing.resize(ntimers, clock_t(0.0));
      Start_time.resize(ntimers, clock_t(0.0));
    }
 
  private:
    Vector<clock_t> Timing;
 
    Vector<clock_t> Start_time;
  };
 
 
  //====================================================================
  //====================================================================
  class DocLinearSolverInfo
  {
  public:
    DocLinearSolverInfo() : Iterations_and_times() {}
 
    void setup_new_time_step()
    {
      // For each new time step, we have a new vector consisting of pairs of
      // unsigned (for the iteration) and double (for the timing results).
      Iterations_and_times.push_back(Vector<Vector<double>>());
    }
 
    void clear_current_time_step()
    {
      Iterations_and_times.pop_back();
      setup_new_time_step();
    }
 
    void add_iteration_and_time(unsigned iter,
                                double prec_setup_time,
                                double linear_solver_time)
    {
#ifdef PARANOID
      if (Iterations_and_times.size() == 0)
      {
        std::ostringstream error_message;
        error_message << "Iterations_and_times is empty. "
                      << "Call setup_new_time_step()\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Vector<double> tmp_iter_time_vec(3, 0);
      tmp_iter_time_vec[0] = (double)iter;
      tmp_iter_time_vec[1] = prec_setup_time;
      tmp_iter_time_vec[2] = linear_solver_time;
      Iterations_and_times.back().push_back(tmp_iter_time_vec);
    }
 
    unsigned current_ntime_step() const
    {
      return Iterations_and_times.size();
    }
 
    unsigned current_nnewton_step() const
    {
      return Iterations_and_times.back().size();
    }
 
    Vector<Vector<Vector<double>>>& iterations_and_times()
    {
      return Iterations_and_times;
    }
 
    Vector<Vector<Vector<double>>> iterations_and_times() const
    {
      return Iterations_and_times;
    }
 
  private:
    Vector<Vector<Vector<double>>> Iterations_and_times;
  };
 
  //====================================================================
  //====================================================================
  class DocInfo
  {
  public:
    DocInfo()
      : Directory("."),
        Doc_flag(true),
        Number(0),
        Label(""),
        Directory_must_exist(true)
    {
    }
 
    DocInfo(const std::string& directory)
      : Directory(directory),
        Doc_flag(true),
        Number(0),
        Label(""),
        Directory_must_exist(true)
    {
    }
 
    std::string directory() const
    {
      return Directory;
    }
 
    void set_directory(const std::string& directory);
 
    void enable_doc()
    {
      Doc_flag = true;
    }
 
    void disable_doc()
    {
      Doc_flag = false;
    }
 
    bool is_doc_enabled() const
    {
      return Doc_flag;
    }
 
    unsigned& number()
    {
      return Number;
    }
 
    unsigned number() const
    {
      return Number;
    }
 
    std::string number_as_string() const
    {
      return StringConversion::to_string(Number);
    }
 
    std::string& label()
    {
      return Label;
    }
 
    std::string label() const
    {
      return Label;
    }
 
    void enable_error_if_directory_does_not_exist()
    {
      Directory_must_exist = true;
    }
 
    void disable_error_if_directory_does_not_exist()
    {
      Directory_must_exist = false;
    }
 
  private:
    std::string Directory;
 
    bool Doc_flag;
 
    unsigned Number;
 
    std::string Label;
 
    bool Directory_must_exist;
  };
 
 
  //====================================================================
  // Namespace for command line arguments
  //====================================================================
  namespace CommandLineArgs
  {
    template<class T>
    struct ArgInfo
    {
      ArgInfo(const bool& is_set, T* arg_pt, const std::string& doc)
      {
        this->is_set = is_set;
        this->arg_pt = arg_pt;
        this->doc = doc;
      }
 
      ArgInfo()
      {
        this->is_set = false;
        this->arg_pt = 0;
        this->doc = "";
      }
 
      // /// Assignment operator. We need this to be able to store these things
      // in
      // /// maps.
      // void operator=(const ArgInfo& that)
      //  {
      //   this->is_set = that.is_set;
      //   this->arg_pt = that.arg_pt;
      //   this->doc = that.doc;
      //  }
 
      bool is_set;
 
      T* arg_pt;
 
      std::string doc;
    };
 
    // Number of arguments + 1
    extern int Argc;
 
    // Arguments themselves
    extern char** Argv;
 
    extern std::map<std::string, ArgInfo<bool>> Specified_command_line_flag;
 
    extern std::map<std::string, ArgInfo<double>>
      Specified_command_line_double_pt;
 
    extern std::map<std::string, ArgInfo<int>> Specified_command_line_int_pt;
 
    extern std::map<std::string, ArgInfo<unsigned>>
      Specified_command_line_unsigned_pt;
 
    extern std::map<std::string, ArgInfo<std::string>>
      Specified_command_line_string_pt;
 
    // Set values
    extern void setup(int argc, char** argv);
 
    // Doc the command line arguments
    extern void output();
 
    extern void specify_command_line_flag(
      const std::string& command_line_flag,
      const std::string& documentation = "Undocumented");
 
    extern void specify_command_line_flag(
      const std::string& command_line_flag,
      double* arg_pt,
      const std::string& documentation = "Undocumented");
 
    extern void specify_command_line_flag(
      const std::string& command_line_flag,
      int* arg_pt,
      const std::string& documentation = "Undocumented");
 
    extern void specify_command_line_flag(
      const std::string& command_line_flag,
      unsigned* arg_pt,
      const std::string& documentation = "Undocumented");
 
    extern void specify_command_line_flag(
      const std::string& command_line_flag,
      std::string* arg_pt,
      const std::string& documentation = "Undocumented");
 
    extern bool command_line_flag_has_been_set(const std::string& flag);
 
    extern void doc_all_flags(
      std::ostream& outstream = *oomph_info.stream_pt());
 
    extern void doc_specified_flags();
 
    extern void doc_available_flags();
 
    extern void check_arg_index(const int& argc, const int& arg_index);
 
    extern void parse_and_assign(
      int argc, char* argv[], const bool& throw_on_unrecognised_args = false);
 
    extern void parse_and_assign(
      const bool& throw_on_unrecognised_args = false);
 
  } // namespace CommandLineArgs
 
  // forward declaration of OomphCommunicator class
  class OomphCommunicator;
 
#ifdef OOMPH_HAS_MPI
  //========================================================================
  //========================================================================
  class MPIOutputModifier : public OutputModifier
  {
  private:
    unsigned Output_rank;
 
    bool Output_from_single_processor;
 
    OomphCommunicator* Communicator_pt;
 
  public:
    MPIOutputModifier() : Output_rank(0), Output_from_single_processor(false) {}
 
    OomphCommunicator*& communicator_pt()
    {
      return Communicator_pt;
    }
 
    virtual bool operator()(std::ostream& stream);
 
    void restrict_output_to_single_processor(const unsigned& output_rank = 0)
    {
      Output_from_single_processor = true;
      Output_rank = output_rank;
    }
 
 
    void allow_output_from_all_processors()
    {
      Output_from_single_processor = false;
    }
  };
 
 
  //========================================================================
  //========================================================================
  extern MPIOutputModifier oomph_mpi_output;
 
  //=== Forward DenseMatrix class
  template<class T>
  class DenseMatrix;
 
  // forward declaration of oomph-communicator class
  // class OomphCommunicator;
 
#endif
 
 
  //======================================================================
  //=====================================================================
  class MPI_Helpers
  {
  public:
    static void init(int argc,
                     char** argv,
                     const bool& make_duplicate_of_mpi_comm_world = true);
 
    static void finalize();
 
    static OomphCommunicator* communicator_pt();
 
    static bool mpi_has_been_initialised()
    {
      return MPI_has_been_initialised;
    }
 
  private:
    MPI_Helpers();
 
    MPI_Helpers(const MPI_Helpers&);
 
    static bool MPI_has_been_initialised;
 
    static OomphCommunicator* Communicator_pt;
  };
 
 
  //====================================================================
  // Namespace for flagging up obsolete parts of the code
  //====================================================================
  namespace ObsoleteCode
  {
    // Flag up obsolete parts of the code
    extern bool FlagObsoleteCode;
 
    // Output warning message
    extern void obsolete();
 
    // Output specified warning message
    extern void obsolete(const std::string& message);
 
  } // namespace ObsoleteCode
 
  //====================================================================
  // Namespace for tecplot stuff
  //====================================================================
  namespace TecplotNames
  {
    // Tecplot colours
    extern Vector<std::string> colour;
 
    // Setup tecplot colours namespace
    extern void setup();
 
  } // namespace TecplotNames
 
 
#ifdef LEAK_CHECK
 
  //====================================================================
  // Namespace for leak check: Keep a running count of all instantiated
  // objects -- add your own if you want to...
  //====================================================================
  namespace LeakCheckNames
  {
    extern long QuadTree_build;
    extern long OcTree_build;
    extern long QuadTreeForest_build;
    extern long OcTreeForest_build;
    extern long RefineableQElement<2> _build;
    extern long RefineableQElement<3> _build;
    extern long MacroElement_build;
    extern long HangInfo_build;
    extern long Node_build;
    extern long GeomReference_build;
    extern long AlgebraicNodeNodeUpdateInfo_build;
    extern long AlgebraicNode_build;
 
    // Reset counters
    extern void reset();
 
    // Doc counters
    extern void doc();
  } // namespace LeakCheckNames
 
#endif
 
  //====================================================================
  // Namespace for pause() command
  //====================================================================
  namespace PauseFlags
  {
    // Flag to enable pausing code
    extern bool PauseFlag;
 
  } // namespace PauseFlags
 
  void pause(std::string message);
 
 
  //=============================================================================
  //=============================================================================
  namespace TimingHelpers
  {
    double timer();
 
    std::string convert_secs_to_formatted_string(const double& time_in_sec);
 
  } // end of namespace TimingHelpers
 
 
 
 
  //===============================================================
  //===============================================================
  namespace MemoryUsage
  {
    extern bool Suppress_mpi_synchronisation;
 
    extern std::string My_memory_usage_system_string;
 
    extern bool Bypass_all_memory_usage_monitoring;
 
    extern std::string My_memory_usage_filename;
 
    void empty_my_memory_usage_file();
 
#ifdef OOMPH_HAS_UNISTDH
 
    void doc_my_memory_usage(const std::string& prefix_string = 0);
 
#endif
 
    extern std::string Total_memory_usage_system_string;
 
    extern std::string Total_memory_usage_filename;
 
    void empty_total_memory_usage_file();
 
    void doc_total_memory_usage(const std::string& prefix_string = "");
 
    void empty_memory_usage_files();
 
    void doc_memory_usage(const std::string& prefix_string = "");
 
    extern std::string Top_system_string;
 
    extern std::string Top_output_filename;
 
    void empty_top_file();
 
    void run_continous_top(const std::string& comment = "");
 
    void stop_continous_top(const std::string& comment = "");
 
    void insert_comment_to_continous_top(const std::string& comment);
 
  } // end of namespace MemoryUsage
 
 
 
  // Forward decl.
  class Problem;
  template<class T>
  class DenseMatrix;
 
  class SolutionFunctorBase
  {
  public:
    // Some typedefs first:
 
    typedef double (*TimeSpaceToDoubleFctPt)(const double& t,
                                             const Vector<double>& x);
 
    typedef Vector<double> (*TimeSpaceToDoubleVectFctPt)(
      const double& t, const Vector<double>& x);
 
    typedef Vector<double> (*TimeSpaceValueToDoubleVectFctPt)(
      const double& t, const Vector<double>& x, const Vector<double>& u);
 
    virtual ~SolutionFunctorBase() {}
 
    virtual Vector<double> operator()(const double& t,
                                      const Vector<double>& x) const = 0;
 
    virtual Vector<double> derivative(const double& t,
                                      const Vector<double>& x,
                                      const Vector<double>& u) const = 0;
 
    virtual void jacobian(const double& t,
                          const Vector<double>& x,
                          const Vector<double>& u,
                          DenseMatrix<double>& jacobian) const
    {
      std::string err = "No Jacobian function implemented";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual bool have_jacobian() const
    {
      return false;
    }
 
    virtual void initialise_from_problem(const Problem* problem_pt) {}
  };
 
 
  class SolutionFunctor : public SolutionFunctorBase
  {
    // This could easily be extended to take
    // FiniteElement::UnsteadyExactSolutionFctPt function (i.e. functions
    // where the output is placed into a given vector rather than returned)
    // as well--just add the appropriate storage pointers and constructors.
 
  public:
    // Constructors:
 
    SolutionFunctor()
    {
      Solution_fpt = 0;
      Derivative_fpt = 0;
    }
 
    SolutionFunctor(TimeSpaceToDoubleVectFctPt solution_fpt)
    {
      Solution_fpt = solution_fpt;
      Derivative_fpt = 0;
    }
 
    SolutionFunctor(TimeSpaceToDoubleVectFctPt solution_fpt,
                    TimeSpaceValueToDoubleVectFctPt derivative_fpt)
    {
      Solution_fpt = solution_fpt;
      Derivative_fpt = derivative_fpt;
    }
 
    virtual ~SolutionFunctor() {}
 
    SolutionFunctor(const SolutionFunctor& that)
    {
      Solution_fpt = that.Solution_fpt;
      Derivative_fpt = that.Derivative_fpt;
    }
 
    void operator=(const SolutionFunctor& that)
    {
      this->Solution_fpt = that.Solution_fpt;
      this->Derivative_fpt = that.Derivative_fpt;
    }
 
    virtual Vector<double> operator()(const double& t,
                                      const Vector<double>& x) const
    {
#ifdef PARANOID
      if (Solution_fpt == 0)
      {
        std::string err = "Solution_fpt is null!";
        throw OomphLibError(
          err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Solution_fpt(t, x);
    }
 
    virtual Vector<double> derivative(const double& t,
                                      const Vector<double>& x,
                                      const Vector<double>& u) const
    {
#ifdef PARANOID
      if (Derivative_fpt == 0)
      {
        std::string err = "Derivative_fpt is null!";
        throw OomphLibError(
          err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Derivative_fpt(t, x, u);
    }
 
    TimeSpaceToDoubleVectFctPt Solution_fpt;
 
    TimeSpaceValueToDoubleVectFctPt Derivative_fpt;
  };
 
 
} // namespace oomph
#endif