Iteration.h

Go to the documentation of this file.

// Copyright (C) 2003 Johan Hoffman and Anders Logg.
// Licensed under the GNU GPL Version 2.
//
// - Updates by Johan Jansson (2004)

#ifndef __ITERATION_H
#define __ITERATION_H

#include <dolfin/dolfin_log.h>
#include <dolfin/constants.h>

namespace dolfin
{
  class Solution;
  class RHS;
  class TimeSlab;
  class Element;
  class ElementGroup;
  class ElementGroupList;
  class FixedPointIteration;


  class Iteration
  {
  public:

    // Type of iteration
    enum State {nonstiff, stiff1, stiff2, stiff3, stiff, newton};

    // Element residuals
    struct Residuals
    {
      Residuals() : r0(0), r1(0), r2(0) {}
      void operator=(real r) { r1 = r2; r2 = r; }
      void reset() { r0 = 0.0; r1 = 0.0; r2 = 0.0; }
      friend LogStream& operator<<(LogStream& stream, const Residuals& r)
      { stream << "residual: " << r.r1 << " --> " << r.r2; return stream; }
      real r0, r1, r2;
    };

    // Increments
    struct Increments
    {
      Increments() : d0(0), d1(0), d2(0), dtot(0) {}
      void operator=(real d) { d0 = d1; d1 = d2; d2 = d; dtot = std::max(dtot, d); }
      void reset() { d0 = 0.0; d1 = 0.0; d2 = 0.0; dtot = 0.0; }
      friend LogStream& operator<<(LogStream& stream, const Increments& d)
      { stream << "increment: " << d.d0 << " --> " << d.d1 << " --> " << d.d2; return stream; }
      real d0, d1, d2, dtot;
    };

    Iteration(Solution& u, RHS& f, FixedPointIteration& fixpoint,
              unsigned int maxiter, real tol, real maxdiv, real maxconv,
              unsigned int depth);
    
    virtual ~Iteration();

    virtual State state() const = 0;

    // Start iteration on group list
    virtual void start(ElementGroupList& list) = 0;

    // Start iteration on element group
    virtual void start(ElementGroup& group) = 0;
    
    // Start iteration on element
    virtual void start(Element& element) = 0;

    virtual void update(ElementGroupList& list, Increments& d) = 0;

    virtual void update(ElementGroup& group, Increments& d) = 0;

    virtual void update(Element& element, Increments& d) = 0;

    virtual void stabilize(ElementGroupList& list, const Increments& d,
                           unsigned int n) = 0;
    
    virtual void stabilize(ElementGroup& group, const Increments& d,
                           unsigned int n) = 0;
    
    virtual void stabilize(Element& element, const Increments& d,
                           unsigned int n) = 0;
    
    virtual bool converged(ElementGroupList& list, const Increments& d, unsigned int n) = 0;

    virtual bool converged(ElementGroup& group, const Increments& d, unsigned int n) = 0;

    virtual bool converged(Element& element, const Increments& d, unsigned int n) = 0;

    virtual bool diverged(ElementGroupList& list, const Increments& d, unsigned int n, State& newstate) = 0;
    
    virtual bool diverged(ElementGroup& group, const Increments& d, unsigned int n, State& newstate) = 0;

    virtual bool diverged(Element& element, const Increments& d, unsigned int n, State& newstate) = 0;

    virtual void report() const = 0;

    void stabilization(real& alpha, unsigned int& m) const;

    unsigned int depth() const;

    void down();

    void up();

    void init(ElementGroup& group);

    // Update initial data for element
    void init(Element& element);

    void reset(ElementGroupList& list);

    void reset(ElementGroup& group);

    // Reset element
    void reset(Element& element);

    // Compute L2 norm of element residual for group list
    real residual(ElementGroupList& list);

    // Compute L2 norm of element residual for element group
    real residual(ElementGroup& group);

    // Compute absolute value of element residual for element
    real residual(Element& element);
    
    // Check if latest iteration was accepted
    bool accept() const;

  protected:

    // An array of values used for Gauss-Jacobi iteration
    struct Values
    {
      Values();
      ~Values();

      void init(unsigned int size);

      real* values;
      unsigned int size;
      unsigned int offset;
    };

    // Stabilization for adaptive iteration
    bool stabilize(const Increments& d, unsigned int n, real r);
    
    // Compute divergence
    real computeDivergence(ElementGroupList& list);

    // Compute divergence
    real computeDivergence(ElementGroup& group);

    // Unstabilized update
    void updateUnstabilized(ElementGroupList& list, Increments& d);

    // Unstabilized update
    void updateUnstabilized(ElementGroup& group, Increments& d);

    // Compute alpha
    real computeAlpha(real rho) const;

    // Compute number of damping steps
    unsigned int computeSteps(real rho) const;
    
    // Initialize data for propagation of initial values
    void initInitialData(real t0);

    // Initialize data to new size
    void initData(Values& values, unsigned int size);

    // Compute size of data
    unsigned int dataSize(ElementGroupList& list);

    // Compute size of data
    unsigned int dataSize(ElementGroup& group);

    // Copy data from group list
    void copyData(ElementGroupList& list, Values& values);

    // Copy data to group list
    void copyData(Values& values, ElementGroupList& list);

    // Copy data from element group
    void copyData(ElementGroup& group, Values& values);
    
    // Copy data to element group
    void copyData(Values& values, ElementGroup& group);

    // Check if the number is between 0 and inf
    bool positive(real number) const;

    //--- Iteration data ---

    Solution& u;
    RHS& f;
    FixedPointIteration& fixpoint;

    unsigned int maxiter;

    real maxdiv;
    real maxconv;
    real tol;

    // Stabilization parameter
    real alpha;

    // Angle of sector, gamma = cos(theta)
    real gamma;

    // Increment at start of stabilizing iterations
    real d0;

    // Number of stabilizing iterations
    unsigned int m;

    // Number of remaining stabilizing iterations
    unsigned int j;

    // Depth of iteration (0,1,2,3)
    unsigned int _depth;

    // Temporary data used to restore element values after iteration
    Values x0;

    // Temporary data used for Gauss-Jacobi iteration
    Values x1;

    // Temporary data used for propagation of initial values
    Values u0;

    // True if we should keep the new values
    bool _accept;

  };

}

#endif

Documentation automatically generated with Doxygen on 10 Sep 2004.