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:
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)
* 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);
/**
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);
*/
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;
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;
};
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)
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();
// 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)
--- /dev/null
+// ---------------------------------------------------------------------
+// $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));
+}
+
+
--- /dev/null
+
+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
projects / dealii-svn.git / commitdiff