diff --git a/BaseLib/DebugTools.h b/BaseLib/DebugTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc86a2d5bebf74df6fbd5fac8fb7da37cff2d87b
--- /dev/null
+++ b/BaseLib/DebugTools.h
@@ -0,0 +1,29 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2014-03-07
+ * \brief Helper tools for debugging
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef DEBUGTOOLS_H
+#define DEBUGTOOLS_H
+
+#include <ostream>
+#include <algorithm>
+#include <vector>
+
+template<typename T>
+std::ostream &operator <<(std::ostream &os, const std::vector<T> &v) {
+ std::copy(v.begin(), v.end(), std::ostream_iterator<T>(os, " "));
+ os << "\n";
+ return os;
+}
+
+#endif //DEBUGTOOLS_H
+
diff --git a/NumLib/CMakeLists.txt b/NumLib/CMakeLists.txt
index f501692df46851771611e02638e4420dc95ef561..dfa76f13248ea5e4e5ca78d4b98a7bbd1e590768 100644
--- a/NumLib/CMakeLists.txt
+++ b/NumLib/CMakeLists.txt
@@ -5,11 +5,16 @@ SET ( SOURCES ${SOURCES_NUMLIB})
GET_SOURCE_FILES(SOURCES_FEM Fem)
SET ( SOURCES ${SOURCES} ${SOURCES_FEM})
+GET_SOURCE_FILES(SOURCES_TIMESTEP TimeStepping)
+GET_SOURCE_FILES(SOURCES_TIMESTEP_ALGORITHMS TimeStepping/Algorithms)
+SET ( SOURCES ${SOURCES} ${SOURCES_TIMESTEP} ${SOURCES_TIMESTEP_ALGORITHMS})
+
INCLUDE_DIRECTORIES (
.
- ../BaseLib
- ../MathLib
- ../MeshLib
+ ${CMAKE_CURRENT_SOURCE_DIR}/../BaseLib
+ ${CMAKE_CURRENT_SOURCE_DIR}/../GeoLib
+ ${CMAKE_CURRENT_SOURCE_DIR}/../MathLib
+ ${CMAKE_CURRENT_SOURCE_DIR}/../MeshLib
)
@@ -22,5 +27,7 @@ TARGET_LINK_LIBRARIES ( NumLib
BaseLib
GeoLib
MathLib
+ MeshLib
+ logog
)
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3dd163733569767075b6c6211745bf8fb5590e9f
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -0,0 +1,62 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "FixedTimeStepping.h"
+
+#include <cmath>
+#include <numeric>
+#include <limits>
+#include <cassert>
+
+namespace NumLib
+{
+
+FixedTimeStepping::FixedTimeStepping(double t0, double tn, const std::vector<double> &vec_all_dt)
+: _t_initial(t0), _t_end(computeEnd(t0, tn, vec_all_dt)), _dt_vector(vec_all_dt), _ts_prev(t0), _ts_current(t0)
+{}
+
+FixedTimeStepping::FixedTimeStepping(double t0, double tn, double dt)
+: _t_initial(t0), _t_end(tn), _dt_vector(static_cast<std::size_t>(std::ceil((tn-t0)/dt)), dt), _ts_prev(t0), _ts_current(t0)
+{}
+
+const TimeStep FixedTimeStepping::getTimeStep() const
+{
+ return _ts_current;
+}
+
+bool FixedTimeStepping::next()
+{
+ // check if last time step
+ if (_ts_current.steps() == _dt_vector.size()
+ || std::abs(_ts_current.current()-_t_end) < std::numeric_limits<double>::epsilon())
+ return false;
+
+ // confirm current time and move to the next if accepted
+ if (accepted())
+ _ts_prev = _ts_current;
+
+ // prepare the next time step info
+ _ts_current = _ts_prev;
+ double dt = _dt_vector[_ts_prev.steps()];
+ if (_ts_prev.current() + dt > _t_end) // upper bound by t_end
+ dt = _t_end - _ts_prev.current();
+ _ts_current += dt;
+
+ return true;
+}
+
+double FixedTimeStepping::computeEnd(double t_initial, double t_end, const std::vector<double> &dt_vector)
+{
+ double t_sum = t_initial + std::accumulate(dt_vector.begin(), dt_vector.end(), 0);
+ return std::min(t_end, t_sum);
+}
+
+} //NumLib
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
new file mode 100644
index 0000000000000000000000000000000000000000..cdf66028faa9c85463c0b90587217ef9cf54c7fb
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
@@ -0,0 +1,99 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef FIXEDTIMESTEPPING_H_
+#define FIXEDTIMESTEPPING_H_
+
+#include <vector>
+
+#include "ITimeStepAlgorithm.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Fixed time stepping algorithm
+ *
+ * This algorithm returns time step size defined by a user priori.
+ */
+class FixedTimeStepping : public ITimeStepAlgorithm
+{
+public:
+
+ /**
+ * Constructor with homogeneous time step size
+ *
+ * A user provides a single time step size \f$\Delta t\f$. Total number of
+ * time steps is calculated by
+ * \f[
+ * n=\frac{t_{\rm end} - t_{\rm initial}}{\Delta t}
+ * \f].
+ *
+ * @param t_initial start time
+ * @param t_end finish time
+ * @param dt uniform time step size
+ */
+ FixedTimeStepping(double t_initial, double t_end, double dt);
+
+ /**
+ * Constructor with user-specified time step sizes
+ *
+ * A user can specify \f$\Delta t\f$ for each time step (i.e. \f$\Delta t_1,
+ * \Delta t_2, ..., \Delta t_n\f$). Time at \f$m\f$ th step is given as
+ * \f[
+ * t_{m}=\sum_{i=1}^m \Delta t_i + t_{\rm initial}
+ * \f].
+ *
+ * @param t_initial start time
+ * @param t_end finish time
+ * @param vec_all_dt a vector of all time steps
+ */
+ FixedTimeStepping(double t_initial, double t_end, const std::vector<double> &vec_all_dt);
+
+ virtual ~FixedTimeStepping() {}
+
+ /// return the beginning of time steps
+ virtual double begin() const {return _t_initial;}
+
+ /// return the end of time steps
+ virtual double end() const {return _t_end;}
+
+ /// return current time step
+ virtual const TimeStep getTimeStep() const;
+
+ /// move to the next time step
+ virtual bool next();
+
+ /// return if current time step is accepted
+ virtual bool accepted() const {return true;}
+
+ /// return a history of time step sizes
+ virtual const std::vector<double>& getTimeStepSizeHistory() const {return _dt_vector; }
+
+private:
+ /// determine true end time
+ static double computeEnd(double t_initial, double t_end, const std::vector<double> &dt_vector);
+
+ /// initial time
+ const double _t_initial;
+ /// end time
+ const double _t_end;
+ /// a vector of time step sizes
+ const std::vector<double> _dt_vector;
+ /// previous time step information
+ TimeStep _ts_prev;
+ /// current time step information
+ TimeStep _ts_current;
+};
+
+} //NumLib
+
+#endif // FIXEDTIMESTEPPING_H_
diff --git a/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..349c0821e0a1afa48fe0a230e066092c473cdbf4
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
@@ -0,0 +1,53 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef ITIMESTEPALGORITHM_H_
+#define ITIMESTEPALGORITHM_H_
+
+#include <vector>
+
+#include "TimeStepping/TimeStep.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Interface of time stepping algorithms
+ *
+ */
+class ITimeStepAlgorithm
+{
+public:
+ /// return the beginning of time steps
+ virtual double begin() const = 0;
+
+ /// return the end of time steps
+ virtual double end() const = 0;
+
+ /// move to the next time step
+ /// \return true if the next step exists
+ virtual bool next() = 0;
+
+ /// return current time step
+ virtual const TimeStep getTimeStep() const = 0;
+
+ /// return if current time step is accepted or not
+ virtual bool accepted() const = 0;
+
+ /// return a history of time step sizes
+ virtual const std::vector<double>& getTimeStepSizeHistory() const = 0;
+
+ virtual ~ITimeStepAlgorithm() {}
+};
+
+} //NumLib
+
+#endif //ITIMESTEPALGORITHM_H_
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5fbddee5b43fca47f3d52622c9f4d9c4147715df
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
@@ -0,0 +1,100 @@
+/**
+ * \author Haibing Shao and Norihiro Watanabe
+ * \date 2013-08-07
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "IterationNumberBasedAdaptiveTimeStepping.h"
+
+#include <limits>
+#include <algorithm>
+#include <cassert>
+
+namespace NumLib
+{
+
+IterationNumberBasedAdaptiveTimeStepping::IterationNumberBasedAdaptiveTimeStepping(double t0, double tn,
+ double min_ts, double max_ts, double initial_ts,
+ const std::vector<std::size_t> &iter_times_vector,
+ const std::vector<double> &multiplier_vector)
+: _t_initial(t0), _t_end(tn), _iter_times_vector(iter_times_vector), _multiplier_vector(multiplier_vector),
+ _min_ts(min_ts), _max_ts(max_ts), _initial_ts(initial_ts), _max_iter(_iter_times_vector.empty() ? 0 : _iter_times_vector.back()),
+ _iter_times(0), _ts_pre(t0), _ts_current(t0), _n_rejected_steps(0)
+{
+ assert(iter_times_vector.size() == multiplier_vector.size());
+}
+
+bool IterationNumberBasedAdaptiveTimeStepping::next()
+{
+ // check current time step
+ if (std::abs(_ts_current.current()-_t_end) < std::numeric_limits<double>::epsilon())
+ return false;
+
+ // confirm current time and move to the next if accepted
+ if (accepted()) {
+ _ts_pre = _ts_current;
+ _dt_vector.push_back(_ts_current.dt());
+ } else {
+ ++_n_rejected_steps;
+ }
+
+ // prepare the next time step info
+ _ts_current = _ts_pre;
+ _ts_current += getNextTimeStepSize();
+
+ return true;
+}
+
+const TimeStep IterationNumberBasedAdaptiveTimeStepping::getTimeStep() const
+{
+ return _ts_current;
+}
+
+double IterationNumberBasedAdaptiveTimeStepping::getNextTimeStepSize() const
+{
+ double dt = 0.0;
+
+ // if this is the first time step
+ // then we use initial guess provided by a user
+ if ( _ts_pre.steps() == 0 )
+ {
+ dt = _initial_ts;
+ }
+ else // not the first time step
+ {
+ double tmp_multiplier = 1.0;
+ // get the first multiplier by default
+ if (!_multiplier_vector.empty())
+ tmp_multiplier = _multiplier_vector[0];
+ // finding the right multiplier
+ for (std::size_t i=0; i<_iter_times_vector.size(); i++ )
+ if ( this->_iter_times > _iter_times_vector[i] )
+ tmp_multiplier = _multiplier_vector[i];
+ // multiply the the multiplier
+ dt = _ts_pre.dt() * tmp_multiplier;
+ }
+
+ // check whether out of the boundary
+ if ( dt < _min_ts )
+ dt = _min_ts;
+ else if ( dt > _max_ts )
+ dt = _max_ts;
+
+ double t_next = dt + _ts_pre.current();
+ if (t_next > end())
+ dt = end() - _ts_pre.current();
+
+ return dt;
+}
+
+bool IterationNumberBasedAdaptiveTimeStepping::accepted() const
+{
+ return ( this->_iter_times <= this->_max_iter );
+}
+
+} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h
new file mode 100644
index 0000000000000000000000000000000000000000..0c297c7ecb490177ab237efa871d9db9f5c3f090
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h
@@ -0,0 +1,145 @@
+/**
+ * \author Haibing Shao and Norihiro Watanabe
+ * \date 2013-08-07
+ *
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.com)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.com/LICENSE.txt
+ */
+
+#ifndef ITERATIONNUMBERBASEDADAPTIVETIMESTEPPING_H_
+#define ITERATIONNUMBERBASEDADAPTIVETIMESTEPPING_H_
+
+#include <vector>
+
+#include "ITimeStepAlgorithm.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Iteration number based adaptive time stepping
+ *
+ * This algorithm estimates a time step size depending on the number of
+ * iterations (e.g. of iterative linear solvers, nonlinear methods, partitioned
+ * coupling) needed in the last time step (see Hoffmann (2010) for Newton-Raphson case).
+ * The new time step \f$\Delta t_{n+1}\f$ size is calculated as
+ * \f[
+ * \Delta t_{n+1} = \alpha \Delta t_n
+ * \f]
+ * with the previous time step size \f$\Delta t_{n}\f$ and a multiplier coefficient
+ * \f$\alpha\f$ depending on the iteration number.
+ * Note that a time step size is always bounded by the minimum and maximum allowed
+ * value.
+ * \f[
+ * \Delta t_{\min} \le \Delta t \le \Delta t_{\max}
+ * \f]
+ *
+ * For example, users can setup the following time stepping strategy based on
+ * the iteration number of the Newton-Raphson method in the last time step.
+ * <table border="1">
+ * <tr><th>Num. of Newton steps</th><th>0-2</th><th>3-6</th><th>7-8</th><th>9<</th></tr>
+ * <tr><th>Time step size multiplier</th><th>1.6</th><th>1.</th><th>0.5</th><th>0.25 (repeat time step)</th></tr>
+ * <tr><th>Upper and lower bound</th><th colspan="4"> \f$ 1. \le \Delta t \le 10.\f$ </th></tr>
+ * </table>
+ * A time step size is increased for the small iteration number, and decreased for the
+ * large iteration number. If the iteration number exceeds a user-defined threshold (e.g. 9),
+ * a time step is repeated with a smaller time step size.
+ *
+ * Reference
+ * - Hoffmann J (2010) Reactive Transport and Mineral Dissolution/Precipitation
+ * in Porous Media:Efficient Solution Algorithms, Benchmark Computations and
+ * Existence of Global Solutions. PhD thesis. pp82. Friedrich-Alexander-Universität Erlangen-Nürnberg.
+ *
+ */
+class IterationNumberBasedAdaptiveTimeStepping : public ITimeStepAlgorithm
+{
+public:
+
+ /**
+ * Constructor
+ *
+ * @param t_initial start time
+ * @param t_end end time
+ * @param min_ts the minimum allowed time step size
+ * @param max_ts the maximum allowed time step size
+ * @param initial_ts initial time step size
+ * @param iter_times_vector a vector of iteration numbers
+ * (\f$i_1\f$, \f$i_2\f$, ..., \f$i_n\f$) which defines intervals as
+ * \f$[i_1,i_2)\f$, \f$[i_2,i_3)\f$, ..., \f$[i_n,\infty)\f$.
+ * If an iteration number is larger than \f$i_n\f$, current time step is
+ * repeated with the new time step size.
+ * @param multiplier_vector a vector of multiplier coefficients
+ * (\f$a_1\f$, \f$a_2\f$, ..., \f$a_n\f$) corresponding to the intervals given by iter_times_vector.
+ * A time step size is calculated by \f$\Delta t_{n+1} = a * \Delta t_{n}\f$
+ */
+ IterationNumberBasedAdaptiveTimeStepping( double t_initial,
+ double t_end,
+ double min_ts,
+ double max_ts,
+ double initial_ts,
+ const std::vector<std::size_t> &iter_times_vector,
+ const std::vector<double> &multiplier_vector);
+
+ virtual ~IterationNumberBasedAdaptiveTimeStepping() {}
+
+ /// return the beginning of time steps
+ virtual double begin() const {return _t_initial;}
+
+ /// return the end of time steps
+ virtual double end() const {return _t_end;}
+
+ /// return current time step
+ virtual const TimeStep getTimeStep() const;
+
+ /// move to the next time step
+ virtual bool next();
+
+ /// return if the current step is accepted
+ virtual bool accepted() const;
+
+ /// return a history of time step sizes
+ virtual const std::vector<double>& getTimeStepSizeHistory() const {return this->_dt_vector;}
+
+ /// set the number of iterations
+ void setNIterations(std::size_t n_itr) {this->_iter_times = n_itr;}
+
+ /// return the number of repeated steps
+ std::size_t getNRepeatedSteps() const {return this->_n_rejected_steps;}
+
+private:
+ /// calculate the next time step size
+ double getNextTimeStepSize() const;
+
+ /// initial time
+ const double _t_initial;
+ /// end time
+ const double _t_end;
+ /// this vector stores the number of iterations to which the respective multiplier coefficient will be applied
+ const std::vector<std::size_t> _iter_times_vector;
+ /// this vector stores the multiplier coefficients
+ const std::vector<double> _multiplier_vector;
+ /// the minimum allowed time step size
+ const double _min_ts;
+ /// the maximum allowed time step size
+ const double _max_ts;
+ /// initial time step size
+ const double _initial_ts;
+ /// the maximum allowed iteration number to accept current time step
+ const std::size_t _max_iter;
+ /// the number of nonlinear iterations
+ std::size_t _iter_times;
+ /// previous time step
+ TimeStep _ts_pre;
+ /// current time step
+ TimeStep _ts_current;
+ /// history of time step sizes
+ std::vector<double> _dt_vector;
+ /// the number of rejected steps
+ std::size_t _n_rejected_steps;
+};
+
+} // NumLib
+
+#endif // NEWTONADAPTIVETIMESTEPPING_H_
diff --git a/NumLib/TimeStepping/TimeStep.h b/NumLib/TimeStepping/TimeStep.h
new file mode 100644
index 0000000000000000000000000000000000000000..4bd91e3290c64b9925c6b9bd47b3059f76a18eda
--- /dev/null
+++ b/NumLib/TimeStepping/TimeStep.h
@@ -0,0 +1,120 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef TIMESTEP_H_
+#define TIMESTEP_H_
+
+#include <cstddef>
+
+namespace NumLib
+{
+
+/**
+ * \brief Time step object
+ *
+ * TimeStep class contains previous time(\f$t_{n}\f$), current time (\f$t_{n+1}\f$),
+ * time step length (\f$\Delta t_{n+1}\f$), and the number of time steps (\f$n+1\f$).
+ */
+class TimeStep
+{
+public:
+ /**
+ * Initialize a time step
+ * @param current_time current time
+ */
+ explicit TimeStep(double current_time)
+ : _previous(current_time), _current(current_time), _dt(_current-_previous), _steps(0) {}
+
+ /**
+ * Initialize a time step
+ * @param previous_time previous time
+ * @param current_time current time
+ * @param n the number of time steps
+ */
+ TimeStep(double previous_time, double current_time, std::size_t n)
+ : _previous(previous_time), _current(current_time), _dt(_current-_previous), _steps(n) {}
+
+ /// copy a time step
+ TimeStep(const TimeStep &src)
+ : _previous(src._previous), _current(src._current), _dt(_current-_previous), _steps(src._steps) {}
+
+ /// copy a time step
+ TimeStep &operator=(const TimeStep &src)
+ {
+ _previous = src._previous;
+ _current = src._current;
+ _dt = src._dt;
+ _steps = src._steps;
+ return *this;
+ }
+
+ /// return a time step incremented by the given time step size
+ TimeStep operator+(const double dt) const
+ {
+ TimeStep t(*this);
+ t += dt;
+ return t;
+ }
+
+ /// increment time step
+ TimeStep &operator+=(const double dt)
+ {
+ _previous = _current;
+ _current += dt;
+ _dt = dt;
+ _steps++;
+ return *this;
+ }
+
+ /// compare current time
+ bool operator<(const TimeStep &t) const {return (_current < t._current);}
+
+ /// compare current time
+ bool operator<(const double &t) const {return (_current < t);}
+
+ /// compare current time
+ bool operator<=(const TimeStep &t) const {return (_current <= t._current);}
+
+ /// compare current time
+ bool operator<=(const double &t) const {return (_current <= t);}
+
+ /// compare current time
+ bool operator==(const TimeStep &t) const {return (_current == t._current);}
+
+ /// compare current time
+ bool operator==(const double &t) const {return (_current == t);}
+
+ /// return previous time step
+ double previous() const {return _previous;}
+
+ /// return current time step
+ double current() const {return _current;}
+
+ /// time step size from _previous
+ double dt() const {return _dt;}
+
+ /// the number of time _steps
+ std::size_t steps() const {return _steps;}
+
+private:
+ /// previous time step
+ double _previous;
+ /// current time step
+ double _current;
+ /// time step size
+ double _dt;
+ /// the number of time steps
+ std::size_t _steps;
+};
+
+} //NumLib
+
+#endif // TIMESTEP_H_
diff --git a/NumLib/dummy.cpp b/NumLib/dummy.cpp
deleted file mode 100644
index 73c6bbe276e41bbf91aa7c690536f47338edf925..0000000000000000000000000000000000000000
--- a/NumLib/dummy.cpp
+++ /dev/null
@@ -1,4 +0,0 @@
-void dummy()
-{
- return;
-}
diff --git a/Tests/NumLib/TestTimeStep.cpp b/Tests/NumLib/TestTimeStep.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cfd2b2978219f93f59a11f37b5c39d7fad8b1172
--- /dev/null
+++ b/Tests/NumLib/TestTimeStep.cpp
@@ -0,0 +1,62 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include <gtest/gtest.h>
+
+#include "NumLib/TimeStepping/TimeStep.h"
+
+TEST(NumLib, TimeStep)
+{
+ // initial
+ NumLib::TimeStep t1(1.);
+ ASSERT_EQ(1., t1.current());
+ ASSERT_EQ(1., t1.previous());
+ ASSERT_EQ(.0, t1.dt());
+ ASSERT_EQ(0u, t1.steps());
+
+ NumLib::TimeStep t2(0., 1., 1);
+ ASSERT_EQ(1., t2.current());
+ ASSERT_EQ(0., t2.previous());
+ ASSERT_EQ(1., t2.dt());
+ ASSERT_EQ(1u, t2.steps());
+
+ // copy
+ NumLib::TimeStep t3(t2);
+ ASSERT_EQ(1., t3.current());
+ ASSERT_EQ(0., t3.previous());
+ ASSERT_EQ(1., t3.dt());
+ ASSERT_EQ(1u, t3.steps());
+
+ NumLib::TimeStep t4 = t2;
+ ASSERT_EQ(1., t4.current());
+ ASSERT_EQ(0., t4.previous());
+ ASSERT_EQ(1., t4.dt());
+ ASSERT_EQ(1u, t4.steps());
+
+ // increment
+ NumLib::TimeStep t5 = t2 + 2.;
+ ASSERT_EQ(3., t5.current());
+ ASSERT_EQ(1., t5.previous());
+ ASSERT_EQ(2., t5.dt());
+ ASSERT_EQ(2u, t5.steps());
+
+ t5 += 3.;
+ ASSERT_EQ(6., t5.current());
+ ASSERT_EQ(3., t5.previous());
+ ASSERT_EQ(3., t5.dt());
+ ASSERT_EQ(3u, t5.steps());
+
+ // comparison
+ ASSERT_TRUE(t2 == t3);
+ ASSERT_TRUE(t5 == 6.);
+ ASSERT_TRUE(t2 < t5);
+ ASSERT_TRUE(t5 < 7.);
+}
diff --git a/Tests/NumLib/TestTimeSteppingFixed.cpp b/Tests/NumLib/TestTimeSteppingFixed.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e77838b74e5688e8a5cd1a49269aedcc8d250c0
--- /dev/null
+++ b/Tests/NumLib/TestTimeSteppingFixed.cpp
@@ -0,0 +1,72 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include <gtest/gtest.h>
+
+#include <vector>
+
+#include "logog/include/logog.hpp"
+
+#include "NumLib/TimeStepping/TimeStep.h"
+#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+
+#include "../TestTools.h"
+#include "TimeSteppingTestingTools.h"
+
+TEST(NumLib, TimeSteppingFixed)
+{
+ // homogeneous dt
+ {
+ NumLib::FixedTimeStepping fixed(1, 31, 10);
+ const std::vector<double> expected_vec_t = {1, 11, 21, 31};
+
+ std::vector<double> vec_t = timeStepping(fixed);
+
+ ASSERT_EQ(expected_vec_t.size(), vec_t.size());
+ ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+ }
+
+ // dt vector (t_end == t0 + sum(dt))
+ {
+ const std::vector<double> fixed_dt = {10, 10, 10};
+ NumLib::FixedTimeStepping fixed(1, 31, 10);
+ const std::vector<double> expected_vec_t = {1, 11, 21, 31};
+
+ std::vector<double> vec_t = timeStepping(fixed);
+
+ ASSERT_EQ(expected_vec_t.size(), vec_t.size());
+ ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+ }
+
+ // dt vector (t_end < t0 + sum(dt))
+ {
+ const std::vector<double> fixed_dt = {5, 10, 20};
+ NumLib::FixedTimeStepping fixed(1, 31, fixed_dt);
+ const std::vector<double> expected_vec_t = {1, 6, 16, 31};
+
+ std::vector<double> vec_t = timeStepping(fixed);
+
+ ASSERT_EQ(expected_vec_t.size(), vec_t.size());
+ ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+ }
+
+ // dt vector (t_end > t0 + sum(dt))
+ {
+ const std::vector<double> fixed_dt = {5, 10, 10};
+ NumLib::FixedTimeStepping fixed(1, 31, fixed_dt);
+ const std::vector<double> expected_vec_t = {1, 6, 16, 26};
+
+ std::vector<double> vec_t = timeStepping(fixed);
+
+ ASSERT_EQ(expected_vec_t.size(), vec_t.size());
+ ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+ }
+}
diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b2927c0f600b60f64c556e12181930102186a78
--- /dev/null
+++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
@@ -0,0 +1,122 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include <gtest/gtest.h>
+
+#include <vector>
+
+#include "logog/include/logog.hpp"
+
+#include "BaseLib/DebugTools.h"
+#include "NumLib/TimeStepping/TimeStep.h"
+#include "NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h"
+
+#include "../TestTools.h"
+#include "TimeSteppingTestingTools.h"
+
+TEST(NumLib, TimeSteppingIterationNumberBased1)
+{
+ std::vector<std::size_t> iter_times_vector = {0, 3, 5, 7};
+ std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
+ NumLib::IterationNumberBasedAdaptiveTimeStepping alg(1, 31, 1, 10, 1, iter_times_vector, multiplier_vector);
+
+ ASSERT_TRUE(alg.next()); // t=2, dt=1
+ NumLib::TimeStep ts = alg.getTimeStep();
+ ASSERT_EQ(1u, ts.steps());
+ ASSERT_EQ(1., ts.previous());
+ ASSERT_EQ(2., ts.current());
+ ASSERT_EQ(1., ts.dt());
+ ASSERT_TRUE(alg.accepted());
+
+ ASSERT_TRUE(alg.next()); // t=4, dt=2
+
+ // dt*=2
+ alg.setNIterations(3);
+ ASSERT_TRUE(alg.next()); // t=8, dt=4
+ ts = alg.getTimeStep();
+ ASSERT_EQ(3u, ts.steps());
+ ASSERT_EQ(4., ts.previous());
+ ASSERT_EQ(8., ts.current());
+ ASSERT_EQ(4., ts.dt());
+ ASSERT_TRUE(alg.accepted());
+
+ // dt*=1
+ alg.setNIterations(5);
+ ASSERT_TRUE(alg.next()); // t=12, dt=4
+ ts = alg.getTimeStep();
+ ASSERT_EQ(4u, ts.steps());
+ ASSERT_EQ(8., ts.previous());
+ ASSERT_EQ(12., ts.current());
+ ASSERT_EQ(4., ts.dt());
+ ASSERT_TRUE(alg.accepted());
+
+ // dt*=0.5
+ alg.setNIterations(7);
+ ASSERT_TRUE(alg.next()); // t=14, dt=2
+ ts = alg.getTimeStep();
+ ASSERT_EQ(5u, ts.steps());
+ ASSERT_EQ(12., ts.previous());
+ ASSERT_EQ(14., ts.current());
+ ASSERT_EQ(2., ts.dt());
+ ASSERT_TRUE(alg.accepted());
+
+ // dt*=0.25 but dt_min = 1
+ alg.setNIterations(8); // exceed max
+ ASSERT_TRUE(alg.next()); // t=13, dt=1
+ ts = alg.getTimeStep();
+ ASSERT_EQ(5u, ts.steps());
+ ASSERT_EQ(12., ts.previous());
+ ASSERT_EQ(13, ts.current());
+ ASSERT_EQ(1., ts.dt());
+ ASSERT_FALSE(alg.accepted());
+
+ // restart, dt*=1
+ alg.setNIterations(4);
+ ASSERT_TRUE(alg.next()); // t=14, dt=1
+ ts = alg.getTimeStep();
+ ASSERT_EQ(6u, ts.steps());
+ ASSERT_EQ(13., ts.previous());
+ ASSERT_EQ(14, ts.current());
+ ASSERT_EQ(1., ts.dt());
+ ASSERT_TRUE(alg.accepted());
+}
+
+TEST(NumLib, TimeSteppingIterationNumberBased2)
+{
+ std::vector<std::size_t> iter_times_vector = {0, 3, 5, 7};
+ std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
+ NumLib::IterationNumberBasedAdaptiveTimeStepping alg(1, 31, 1, 10, 1, iter_times_vector, multiplier_vector);
+
+ std::vector<std::size_t> nr_iterations = {2, 2, 2, 4, 6, 8, 4, 4, 2, 2};
+ const std::vector<double> expected_vec_t = {1, 2, 4, 8, 16, 24, 26, 28, 30, 31};
+
+ struct IterationNumberUpdate
+ {
+ IterationNumberUpdate(const std::vector<std::size_t> &vec) : _nr_iterations(vec) {}
+ std::vector<std::size_t> _nr_iterations;
+ void operator()(NumLib::IterationNumberBasedAdaptiveTimeStepping &obj)
+ {
+ static std::size_t i = 0;
+ std::size_t n = (i<_nr_iterations.size()) ? _nr_iterations[i++] : 0;
+ //INFO("-> NR-iterations=%d", n);
+ obj.setNIterations(n);
+ }
+ };
+
+ IterationNumberUpdate update(nr_iterations);
+
+ std::vector<double> vec_t = timeStepping(alg, &update);
+ //std::cout << vec_t;
+
+ ASSERT_EQ(expected_vec_t.size(), vec_t.size());
+ ASSERT_EQ(1u, alg.getNRepeatedSteps());
+ ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+}
diff --git a/Tests/NumLib/TimeSteppingTestingTools.h b/Tests/NumLib/TimeSteppingTestingTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f5d11e9f49bfe98a8fe1fdc0c5419c2a39f600f
--- /dev/null
+++ b/Tests/NumLib/TimeSteppingTestingTools.h
@@ -0,0 +1,51 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef TIMESTEPPINGTESTINGTOOLS_H_
+#define TIMESTEPPINGTESTINGTOOLS_H_
+
+#include "logog/include/logog.hpp"
+
+#include "NumLib/TimeStepping/TimeStep.h"
+
+
+namespace
+{
+
+struct Dummy
+{
+ template <class T>
+ void operator()(T &/*obj*/) {}
+};
+
+template <class T_TIME_STEPPING, class T=Dummy>
+std::vector<double> timeStepping(T_TIME_STEPPING &algorithm, T* obj=nullptr)
+{
+ std::vector<double> vec_t;
+ vec_t.push_back(algorithm.begin());
+
+ while (algorithm.next()) {
+ NumLib::TimeStep t = algorithm.getTimeStep();
+ //INFO("t: n=%d,t=%g,dt=%g", t.steps(), t.current(), t.dt());
+ if (obj)
+ (*obj)(algorithm); // do something
+ if (algorithm.accepted()) {
+ vec_t.push_back(t.current());
+ } else {
+ //INFO("*** rejected.");
+ }
+ }
+
+ return vec_t;
+}
+} // namespace
+
+#endif // TIMESTEPPINGTESTINGTOOLS_H_