implicit_midpoint_rule.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//====================================================================
#ifndef OOMPH_IMPLICIT_MIDPOINT_RULE_H
#define OOMPH_IMPLICIT_MIDPOINT_RULE_H
 
// oomph-lib headers
#include "Vector.h"
#include "nodes.h"
#include "matrices.h"
#include "timesteppers.h"
#include "explicit_timesteppers.h"
 
namespace oomph
{
  // Forward decl. so that we can have function of Problem*
  class Problem;
 
 
  class IMRBase : public TimeStepper
  {
  public:
    IMRBase(const bool& adaptive = false)
      : TimeStepper(2, 1) // initialise weights later
    {
      Adaptive_Flag = adaptive;
      Is_steady = false;
      Type = "Midpoint method";
      Predict_by_explicit_step = true;
 
      // If adaptive then we are storing predicted values in slot
      // 4. Otherwise we aren't storing them so leave it as -1.
      if (adaptive)
      {
        Predictor_storage_index = 4;
      }
 
 
      // Storage for weights needs to be 2x(2 + 0/2) (the +2 is extra history
      // for ebdf3 predictor if adaptive). This means we provide ways to
      // calculate the zeroth and first derivatives and in calculations we
      // use 2 + 0/2 time values.
 
      // But actually (for some reason...) the amount of data storage
      // allocated for the timestepper in each node is determined by the
      // size of weight. So we need to store an extra dummy entry in order
      // to have storage space for the predicted value at t_{n+1}.
 
      // Just leave space for predictor values anyway, it's not expensive.
      Weight.resize(2, 5, 0.0);
 
      // Assume that set_weights() or make_steady() is called before use to
      // set up the values stored in Weight.
    }
 
    virtual ~IMRBase() {}
 
    virtual void set_weights() = 0;
 
    virtual unsigned nprev_values_for_value_at_evaluation_time() const = 0;
 
    unsigned order() const
    {
      return 2;
    }
 
    unsigned ndt() const
    {
      return nprev_values();
    }
 
    unsigned nprev_values() const
    {
      return 4;
    }
 
    void shift_time_values(Data* const& data_pt);
 
    void shift_time_positions(Node* const& node_pt);
 
    void set_error_weights() {}
 
    void set_predictor_weights() {}
 
    void assign_initial_values_impulsive(Data* const& data_pt)
    {
      std::string err = "Not implemented";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    void assign_initial_positions_impulsive(Node* const& node_pt)
    {
      std::string err = "Not implemented";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void calculate_predicted_positions(Node* const& node_pt)
    {
      std::string err = "Not implemented";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    double temporal_error_in_position(Node* const& node_pt, const unsigned& i)
    {
      std::string err = "Not implemented";
      throw OomphLibError(
        err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Adaptivity
    void calculate_predicted_values(Data* const& data_pt);
    double temporal_error_in_value(Data* const& data_pt, const unsigned& i);
  };
 
  class IMR : public IMRBase
  {
  public:
 
 
    IMR(const bool& adaptive = false) : IMRBase(adaptive) {}
 
    virtual ~IMR() {}
 
    void set_weights()
    {
      // Set the weights for zero-th derivative (i.e. the value to use in
      // newton solver calculations, implicit midpoint method uses the
      // average of previous and current values).
      Weight(0, 0) = 0.5;
      Weight(0, 1) = 0.5;
 
      // Set the weights for 1st time derivative
      double dt = Time_pt->dt(0);
      Weight(1, 0) = 1.0 / dt;
      Weight(1, 1) = -1.0 / dt;
    }
 
    unsigned nprev_values_for_value_at_evaluation_time() const
    {
      return 2;
    }
 
  private:
    IMR(const IMR& dummy) {}
 
    void operator=(const IMR& dummy) {}
  };
 
 
  class IMRByBDF : public IMRBase
  {
  public:
    IMRByBDF(const bool& adaptive = false) : IMRBase(adaptive)
    {
      Update_pinned = true;
    }
 
    virtual ~IMRByBDF() {}
 
    void set_weights()
    {
      // Use weights from bdf1
      double dt = Time_pt->dt(0);
      Weight(1, 0) = 1.0 / dt;
      Weight(1, 1) = -1.0 / dt;
    }
 
    unsigned nprev_values_for_value_at_evaluation_time() const
    {
      return 1;
    }
 
    void actions_before_timestep(Problem* problem_pt);
 
    void actions_after_timestep(Problem* problem_pt);
 
    bool Update_pinned;
 
  private:
    IMRByBDF(const IMRByBDF& dummy) {}
 
    void operator=(const IMRByBDF& dummy) {}
  };
 
} // namespace oomph
 
#endif