From 6913e1f55ca6c96d0109070d30b76c64274b1676 Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 10 Mar 2014 16:06:58 +0000
Subject: [PATCH] adjust ThetaTimestepping

git-svn-id: https://svn.dealii.org/trunk@32639 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/algorithms/operator.h |   7 +-
 .../deal.II/algorithms/operator.templates.h   |  36 ++++
 .../deal.II/algorithms/theta_timestepping.h   | 117 +++++-------
 .../algorithms/theta_timestepping.templates.h |  31 ++--
 tests/algorithms/theta_01_compat.cc           | 166 ++++++++++++++++++
 tests/algorithms/theta_01_compat.output       |  13 ++
 6 files changed, 287 insertions(+), 83 deletions(-)
 create mode 100644 tests/algorithms/theta_01_compat.cc
 create mode 100644 tests/algorithms/theta_01_compat.output

diff --git a/deal.II/include/deal.II/algorithms/operator.h b/deal.II/include/deal.II/algorithms/operator.h
index 0bc6575e77..6254d35019 100644
--- a/deal.II/include/deal.II/algorithms/operator.h
+++ b/deal.II/include/deal.II/algorithms/operator.h
@@ -196,7 +196,12 @@ namespace Algorithms
     OutputOperator<VECTOR> &operator<< (unsigned int step);
 
     /**
-     * Output all the vectors in NamedData.
+     * Output all the vectors in AnyData.
+     */
+    virtual OutputOperator<VECTOR> &operator<< (const AnyData& vectors);
+    
+    /**
+     * @deprecated Output all the vectors in NamedData.
      */
     virtual OutputOperator<VECTOR> &operator<< (const NamedData<VECTOR *> &vectors);
   protected:
diff --git a/deal.II/include/deal.II/algorithms/operator.templates.h b/deal.II/include/deal.II/algorithms/operator.templates.h
index e2a7458569..5b16380492 100644
--- a/deal.II/include/deal.II/algorithms/operator.templates.h
+++ b/deal.II/include/deal.II/algorithms/operator.templates.h
@@ -119,6 +119,42 @@ namespace Algorithms
     os =&stream;
   }
 
+  template <class VECTOR>
+  OutputOperator<VECTOR> &
+  OutputOperator<VECTOR>::operator<< (const AnyData& vectors)
+  {
+    if (os == 0)
+      {
+        //TODO: make this possible
+        //deallog << ' ' << step;
+        //for (unsigned int i=0;i<vectors.size();++i)
+        //  vectors(i)->print(deallog);
+        //deallog << std::endl;
+      }
+    else
+      {
+        (*os) << ' ' << step;
+        for (unsigned int i=0; i<vectors.size(); ++i)
+	  {
+	    if (vectors.is_type<VECTOR*>(i))
+	      {
+		const VECTOR& v = *vectors.entry<VECTOR*>(i);
+		for (unsigned int j=0; j<v.size(); ++j)
+		  (*os) << ' ' << v(j);
+	      }
+	    else if (vectors.is_type<const VECTOR*>(i))
+	      {
+		const VECTOR& v = *vectors.entry<const VECTOR*>(i);
+		for (unsigned int j=0; j<v.size(); ++j)
+		  (*os) << ' ' << v(j);
+	      }
+	  }
+        (*os) << std::endl;
+      }
+    return *this;
+  }
+  
+
   template <class VECTOR>
   OutputOperator<VECTOR> &
   OutputOperator<VECTOR>::operator<< (const NamedData<VECTOR *> &vectors)
diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.h b/deal.II/include/deal.II/algorithms/theta_timestepping.h
index e1a2d0c40d..43ade372df 100644
--- a/deal.II/include/deal.II/algorithms/theta_timestepping.h
+++ b/deal.II/include/deal.II/algorithms/theta_timestepping.h
@@ -188,6 +188,11 @@ namespace Algorithms
      * should contain the initial
      * value in first position. <tt>out</tt>
      */
+    virtual void operator() (AnyData &out, const AnyData &in);
+
+    /**
+     * @deprecated Use  the function with AnyData arguments
+     */
     virtual void operator() (NamedData<VECTOR *> &out, const NamedData<VECTOR *> &in);
 
     /**
@@ -197,11 +202,8 @@ namespace Algorithms
     virtual void notify(const Event &);
 
     /**
-     * Define an operator which
-     * will output the result in
-     * each step. Note that no
-     * output will be generated
-     * without this.
+     * Define an operator which will output the result in each
+     * step. Note that no output will be generated without this.
      */
     void set_output(OutputOperator<VECTOR> &output);
 
@@ -217,34 +219,23 @@ namespace Algorithms
      */
     const double &step_size() const;
     /**
-     * The weight between implicit
-     * and explicit part.
+     * The weight between implicit and explicit part.
      */
     const double &theta() const;
 
     /**
-     * The data handed to the
-     * #op_explicit time stepping
-     * operator.
+     * The data handed to the #op_explicit time stepping operator.
      *
-     * The time in here is the time
-     * at the beginning of the
-     * current step, the time step
-     * is (1-#theta) times the
-     * actual time step.
+     * The time in here is the time at the beginning of the current
+     * step, the time step is (1-#theta) times the actual time step.
      */
     const TimestepData &explicit_data() const;
 
     /**
-     * The data handed to the
-     * #op_implicit time stepping
-     * operator.
+     * The data handed to the #op_implicit time stepping operator.
      *
-     * The time in here is the time
-     * at the beginning of the
-     * current step, the time step
-     * is #theta times the
-     * actual time step.
+     * The time in here is the time at the beginning of the current
+     * step, the time step is #theta times the actual time step.
      */
     const TimestepData &implicit_data() const;
 
@@ -255,84 +246,61 @@ namespace Algorithms
 
   private:
     /**
-     * The object controlling the
-     * time step size and computing
-     * the new time in each step.
+     * The object controlling the time step size and computing the new
+     * time in each step.
      */
     TimestepControl control;
 
     /**
-     * The control parameter theta in the
-     * range <tt>[0,1]</tt>.
+     * The control parameter theta in the range <tt>[0,1]</tt>. It
+     * defaults to 0.5.
      */
     double vtheta;
     /**
-     * Use adaptive #theta if
-     * <tt>true</tt>.
+     * Use adaptive #theta if <tt>true</tt>. Not yet implemented.
      */
     bool adaptive;
 
     /**
-     * The data for the explicit
-     * part of the scheme.
+     * The data for the explicit part of the scheme.
      */
     TimestepData d_explicit;
 
     /**
-     * The data for the implicit
-     * part of the scheme.
+     * The data for the implicit part of the scheme.
      */
     TimestepData d_implicit;
 
 
     /**
-     * The operator computing the
-     * explicit part of the
-     * scheme. This will receive in
-     * its input data the value at
-     * the current time with name
-     * "Current time solution". It
-     * should obtain the current
-     * time and time step size from
-     * explicit_data().
+     * The operator computing the explicit part of the scheme. This
+     * will receive in its input data the value at the current time
+     * with name "Current time solution". It should obtain the current
+     * time and time step size from explicit_data().
      *
-     * Its return value is
-     * <i>Mu+cAu</i>, where
-     * <i>u</i> is the current
-     * state vector, <i>M</i> the
-     * mass matrix, <i>A</i> the
-     * operator in space and
-     * <i>c</i> is the time step
-     * size in explicit_data().
+     * Its return value is <i>Mu+cAu</i>, where <i>u</i> is the
+     * current state vector, <i>M</i> the mass matrix, <i>A</i> the
+     * operator in space and <i>c</i> is the time step size in
+     * explicit_data().
      */
     SmartPointer<Operator<VECTOR>, ThetaTimestepping<VECTOR> > op_explicit;
 
     /**
-     * The operator solving the
-     * implicit part of the
-     * scheme. It will receive in
-     * its input data the vector
-     * "Previous time". Information on the
-     * timestep should be obtained
-     * from implicit_data().
+     * The operator solving the implicit part of the scheme. It will
+     * receive in its input data the vector "Previous
+     * time". Information on the timestep should be obtained from
+     * implicit_data().
      *
-     * Its return value is the
-     * solution <i>u</i> of
-     * <i>Mu-cAu=f</i>, where
-     * <i>f</i> is the dual space
-     * vector found in the
-     * "Previous time" entry of the
-     * input data, <i>M</i> the
-     * mass matrix, <i>A</i> the
-     * operator in space and
-     * <i>c</i> is the time step
+     * Its return value is the solution <i>u</i> of <i>Mu-cAu=f</i>,
+     * where <i>f</i> is the dual space vector found in the "Previous
+     * time" entry of the input data, <i>M</i> the mass matrix,
+     * <i>A</i> the operator in space and <i>c</i> is the time step
      * size in explicit_data().
      */
     SmartPointer<Operator<VECTOR>, ThetaTimestepping<VECTOR> > op_implicit;
 
     /**
-     * The operator writing the
-     * output in each time step
+     * The operator writing the output in each time step
      */
     SmartPointer<OutputOperator<VECTOR>, ThetaTimestepping<VECTOR> > output;
   };
@@ -355,6 +323,15 @@ namespace Algorithms
     return d_implicit;
   }
 
+
+  template <class VECTOR>
+  inline
+  TimestepControl &
+  ThetaTimestepping<VECTOR>::timestep_control ()
+  {
+    return control;
+  }
+
   template <class VECTOR>
   inline
   void ThetaTimestepping<VECTOR>::set_output (OutputOperator<VECTOR> &out)
diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
index b0f1366641..6194088cfe 100644
--- a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
+++ b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
@@ -64,13 +64,23 @@ namespace Algorithms
   template <class VECTOR>
   void
   ThetaTimestepping<VECTOR>::operator() (NamedData<VECTOR *> &out, const NamedData<VECTOR *> &in)
+  {
+    Operator<VECTOR>::operator() (out, in);
+  }
+  
+  
+  template <class VECTOR>
+  void
+  ThetaTimestepping<VECTOR>::operator() (AnyData &out, const AnyData &in)
   {
     Assert(!adaptive, ExcNotImplemented());
 
     deallog.push ("Theta");
+
+    VECTOR& solution = *out.entry<VECTOR*>(0);
     GrowingVectorMemory<VECTOR> mem;
     typename VectorMemory<VECTOR>::Pointer aux(mem);
-    aux->reinit(*out(0));
+    aux->reinit(solution);
 
     control.restart();
 
@@ -79,20 +89,17 @@ namespace Algorithms
     // The data used to compute the
     // vector associated with the old
     // timestep
-    VECTOR *p = out(0);
-    NamedData<VECTOR *> src1;
-    src1.add(p, "Previous iterate");
+    AnyData src1;
+    src1.add<const VECTOR*>(&solution, "Previous iterate");
     src1.merge(in);
 
-    NamedData<VECTOR *> src2;
-    src2.add(p, "Previous iterate");
+    AnyData src2;
+    src2.add<const VECTOR*>(&solution, "Previous iterate");
 
-    p = aux;
-    NamedData<VECTOR *> out1;
-    out1.add(p, "Result");
-    // The data provided to the inner
-    // solver
-    src2.add(p, "Previous time");
+    AnyData out1;
+    out1.add<VECTOR*>(aux, "Solution");
+    // The data provided to the inner solver
+    src2.add<const VECTOR*>(aux, "Previous time");
     src2.merge(in);
 
     if (output != 0)
diff --git a/tests/algorithms/theta_01_compat.cc b/tests/algorithms/theta_01_compat.cc
new file mode 100644
index 0000000000..84041252a2
--- /dev/null
+++ b/tests/algorithms/theta_01_compat.cc
@@ -0,0 +1,166 @@
+// ---------------------------------------------------------------------
+// $Id: theta_timestepping.cc 30050 2013-07-18 20:03:31Z maier $
+//
+// Copyright (C) 2005 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// See documentation of ThetaTimestepping for documentation of this example
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/algorithms/operator.h>
+#include <deal.II/algorithms/theta_timestepping.h>
+
+#include <iostream>
+
+using namespace dealii;
+using namespace Algorithms;
+
+
+class Explicit
+  : public Operator<Vector<double> >
+{
+public:
+  Explicit(const FullMatrix<double> &matrix);
+  void operator() (NamedData<Vector<double>*> &out,
+                   const NamedData<Vector<double>*> &in);
+
+  void initialize_timestep_data(const TimestepData &);
+private:
+  const TimestepData *timestep_data;
+  SmartPointer<const FullMatrix<double>, Explicit> matrix;
+  FullMatrix<double> m;
+};
+
+
+class Implicit
+  : public Operator<Vector<double> >
+{
+public:
+  Implicit(const FullMatrix<double> &matrix);
+  void operator() (NamedData<Vector<double>*> &out,
+                   const NamedData<Vector<double>*> &in);
+
+  void initialize_timestep_data(const TimestepData &);
+private:
+  const TimestepData *timestep_data;
+  SmartPointer<const FullMatrix<double>, Implicit> matrix;
+  FullMatrix<double> m;
+};
+
+// End of declarations
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  FullMatrix<double> matrix(2);
+  matrix(0,0) = 0.;
+  matrix(1,1) = 0.;
+  matrix(0,1) = numbers::PI;
+  matrix(1,0) = -numbers::PI;
+
+  OutputOperator<Vector<double> > out;
+  out.initialize_stream(logfile);
+
+  Explicit op_explicit(matrix);
+  Implicit op_implicit(matrix);
+  ThetaTimestepping<Vector<double> > solver(op_explicit, op_implicit);
+  op_explicit.initialize_timestep_data(solver.explicit_data());
+  op_implicit.initialize_timestep_data(solver.implicit_data());
+  solver.timestep_control().start_step(0.1);
+  //solver.set_output(out);
+
+  Vector<double> value(2);
+  value(0) = 1.;
+  NamedData<Vector<double>*> indata;
+  NamedData<Vector<double>*> outdata;
+  Vector<double> *p = &value;
+  outdata.add(p, "value");
+  deallog << "Initial: " << value(0) << ' ' << value(1) << std::endl;
+  solver.notify(Events::initial);
+  solver(outdata, indata);
+  deallog << "Result: " << value(0) << ' ' << value(1)
+	  << " Norm " << value.l2_norm() << std::endl;
+}
+
+
+Explicit::Explicit(const FullMatrix<double> &M)
+  :
+  matrix(&M)
+{
+  m.reinit(M.m(), M.n());
+}
+
+
+void
+Explicit::initialize_timestep_data(const TimestepData &t)
+{
+  timestep_data = &t;
+}
+
+
+void
+Explicit::operator() (NamedData<Vector<double>*> &out, const NamedData<Vector<double>*> &in)
+{
+  if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+    {
+      m.equ(-timestep_data->step, *matrix);
+      for (unsigned int i=0; i<m.m(); ++i)
+        m(i,i) += 1.;
+    }
+  this->notifications.clear();
+  unsigned int i = in.find("Previous iterate");
+  m.vmult(*out(0), *in(i));
+}
+
+
+Implicit::Implicit(const FullMatrix<double> &M)
+  :
+  matrix(&M)
+{
+  m.reinit(M.m(), M.n());
+}
+
+
+void
+Implicit::initialize_timestep_data(const TimestepData &t)
+{
+  timestep_data = &t;
+}
+
+
+void
+Implicit::operator() (NamedData<Vector<double>*> &out, const NamedData<Vector<double>*> &in)
+{
+  if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+    {
+      m.equ(timestep_data->step, *matrix);
+      for (unsigned int i=0; i<m.m(); ++i)
+        m(i,i) += 1.;
+      m.gauss_jordan();
+    }
+  this->notifications.clear();
+
+  unsigned int i = in.find("Previous time");
+  m.vmult(*out(0), *in(i));
+}
+
+
diff --git a/tests/algorithms/theta_01_compat.output b/tests/algorithms/theta_01_compat.output
new file mode 100644
index 0000000000..7c4d9ddfc7
--- /dev/null
+++ b/tests/algorithms/theta_01_compat.output
@@ -0,0 +1,13 @@
+
+DEAL::Initial: 1.00000 0
+DEAL:Theta::Time step:0.100000
+DEAL:Theta::Time step:0.200000
+DEAL:Theta::Time step:0.300000
+DEAL:Theta::Time step:0.400000
+DEAL:Theta::Time step:0.500000
+DEAL:Theta::Time step:0.600000
+DEAL:Theta::Time step:0.700000
+DEAL:Theta::Time step:0.800000
+DEAL:Theta::Time step:0.900000
+DEAL:Theta::Time step:1.00000
+DEAL::Result: -0.999676 0.0254599 Norm 1.00000
-- 
2.39.5