]> https://gitweb.dealii.org/ - dealii.git/commitdiff
start algorithms library
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 12 Jan 2010 16:56:38 +0000 (16:56 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 12 Jan 2010 16:56:38 +0000 (16:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@20350 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
deal.II/deal.II/include/numerics/event.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker.h
deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/include/numerics/mesh_worker_info.h
deal.II/deal.II/include/numerics/mesh_worker_loop.h
deal.II/deal.II/include/numerics/newton.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/newton.templates.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/operator.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/operator.templates.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/event.all_dimensions.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker.cc
deal.II/deal.II/source/numerics/operators.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/operators.inst.in [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/numerics/event.h b/deal.II/deal.II/include/numerics/event.h
new file mode 100644 (file)
index 0000000..7074e00
--- /dev/null
@@ -0,0 +1,287 @@
+//---------------------------------------------------------------------------
+//    $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__event_h
+#define __deal2__event_h
+
+#include <base/config.h>
+
+#include <vector>
+#include <string>
+#include <iostream>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+/**
+ * Objects of this kind are used to notify interior applications of
+ * changes provoked by an outer loop. They are handed to the
+ * application through Operator::notify() and it is up to the
+ * actual application how to handle them.
+ *
+ * Event is organized as an extensible binary enumerator. Every class
+ * can add its own events using assign(). A typical code example is
+ *
+ * <pre>
+ * class A
+ * {
+ *   static Event event;
+ * };
+ *
+ * Event A::event = Event::assign("Event for A");
+ * </pre>
+ */
+  class Event
+  {
+    public:
+                                      /**
+                                       * This function registers a
+                                       * new event type and assigns a
+                                       * unique identifier to it. The
+                                       * result of this function
+                                       * should be stored for later
+                                       * use.
+                                       */
+      static Event assign (const char* name);
+
+                                      /**
+                                       * If you forgot to store the
+                                       * result of assign, here is
+                                       * how to retrieve it knowing
+                                       * the name.
+                                       */
+//      static Event find(const std::string& name);
+      
+                                      /**
+                                       * Constructor, generating a
+                                       * clear Event.
+                                       */
+      Event ();
+
+                                      /**
+                                       * Clear all flags
+                                       */
+      void clear();
+
+                                      /**
+                                       * Set all flags
+                                       */
+      void all();
+
+                                      /**
+                                       * Add the flags of the other event
+                                       */
+      Event& operator += (const Event& event);
+
+                                      /**
+                                       * Clear the flags of the other event
+                                       */
+      Event& operator -= (const Event& event);
+
+                                      /**
+                                       * Test whether all the flags
+                                       * set in the other Event are
+                                       * also set in this one.
+                                       */
+      bool test (const Event& event) const;
+
+                                      /**
+                                       * Return <tt>true</tt> if any
+                                       * event is set.
+                                       */
+      bool any () const;
+
+                                      /**
+                                       * List the flags to a stream.
+                                       */
+      template <class OS>
+      void print (OS& os) const;
+
+                                      /**
+                                       * List all assigned events.
+                                       */
+      template <class OS>
+      static void print_assigned (OS& os);
+      
+    private:
+                                      /**
+                                       * Sometimes, actions have to
+                                       * be taken by all
+                                       * means. Therefore, if this
+                                       * value is true, test() always
+                                       * returns true.
+                                       */
+      bool all_true;
+      
+                                      /**
+                                       * The actual list of events
+                                       */
+      std::vector<bool> flags;
+      
+                                      /**
+                                       * The names of registered events
+                                       */
+      static std::vector<std::string> names;
+  };
+
+/// Events used by library operators  
+  namespace Events
+  {
+/// The current derivative leads to slow convergence of Newton's method
+    extern const Event bad_derivative;
+/// The time stepping scheme starts a new time step
+    extern const Event new_time;
+/// The time stepping scheme has changed the time step size
+    extern const Event new_time_step_size;
+  }
+  
+
+//----------------------------------------------------------------------//
+
+
+  inline
+  bool
+  Event::any () const
+  {
+    if (all_true) return true;
+    for (std::vector<bool>::const_iterator i=flags.begin();
+        i != flags.end(); ++i)
+      if (*i) return true;
+    return false;
+  }
+
+
+  inline
+  bool
+  Event::test (const Event& event) const
+  {
+    
+                                    // First, test all_true in this
+    if (all_true) return true;
+    
+    const unsigned int n = flags.size();
+    const unsigned int m = event.flags.size();
+    const unsigned int n_min = std::min(n, m);
+    
+                                    // Now, if all_true set in the
+                                    // other, then all must be true
+                                    // in this
+    if (event.all_true)
+      {
+                                        // Non existing flags are
+                                        // always assumed false
+       if (m > n)
+         return false;
+
+                                        // Test all flags separately
+                                        // and return false if one is
+                                        // not set
+       for (std::vector<bool>::const_iterator i=flags.begin();
+            i != flags.end(); ++i)
+         if (!*i) return false;
+                                        // All flags are set
+       return true;
+      }
+
+                                    // Finally, compare each flag
+                                    // separately
+    for (unsigned int i=0;i<n_min;++i)
+      if (event.flags[i] && !flags[i])
+       return false;
+    for (unsigned int i=n_min;i<m;++i)
+      if (event.flags[i])
+       return false;
+    return true;
+  }
+  
+
+  
+  inline
+  Event& Event::operator += (const Event& event)
+  {
+    all_true |= event.all_true;
+    if (all_true) return *this;
+    
+    if (flags.size() < event.flags.size())
+      flags.resize(event.flags.size());
+    for (unsigned int i=0;i<flags.size();++i)
+      flags[i] = flags[i] || event.flags[i];
+
+    return *this;
+  }
+  
+
+  inline
+  Event& Event::operator -= (const Event& event)
+  {
+    if (!event.any()) return *this;
+    
+    all_true = false;
+    if(event.all_true)
+      {
+       for (std::vector<bool>::iterator i=flags.begin();
+            i != flags.end(); ++i)
+         *i = false;
+       return *this;
+      }
+    
+    if (flags.size() < event.flags.size())
+      flags.resize(event.flags.size());
+    for (unsigned int i=0;i<flags.size();++i)
+      if (event.flags[i]) flags[i] = false;
+    
+    return *this;
+  }
+  
+
+  template <class OS>
+  inline
+  void
+  Event::print (OS& os) const
+  {
+    if (all_true)
+      os << " ALL";
+    
+    for (unsigned int i=0;i<flags.size();++i)
+      if (flags[i])
+       os << ' ' << names[i];
+  }
+  
+  
+  template <class OS>
+  inline
+  void
+  Event::print_assigned (OS& os)
+  {
+    for (unsigned int i=0;i<names.size();++i)
+      os << i << '\t' << names[i] << std::endl;
+  }
+  
+  
+                                  /**
+                                   * Output shift operator for
+                                   * events. Calls Event::print().
+                                   *
+                                   * @relates Event
+                                   */
+  template <class OS>
+  OS& operator << (OS& o, const Event& e)
+  {
+    e.print(o);
+    return o;
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 5f956c2f3fe70085412eb1e1d04f3dd9a214a29e..535a2216a93c1202fa953b7a7a2c5a69044bf922 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index c822c195bbef535315eccdd56a371d8319bb46a7..f330e125c7a98646acc5b6e0b5b5551306694cc4 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index 795410beceb21890ff6ad2814cfc1aaf750431fc..853c7496e2f38536c093c3be2e0bee0b8182c222 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index 2d79b8f10a66222d084f3d8f84ec7b29b0b9a290..909af2f9bd64b0b5b62126ea34562dd4bb0f830b 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
diff --git a/deal.II/deal.II/include/numerics/newton.h b/deal.II/deal.II/include/numerics/newton.h
new file mode 100644 (file)
index 0000000..5bb44de
--- /dev/null
@@ -0,0 +1,161 @@
+//---------------------------------------------------------------------------
+//    $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__newton_step_control_h
+#define __deal2__newton_step_control_h
+
+#include <base/smartpointer.h>
+#include <lac/solver_control.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+class ParameterHandler;
+
+namespace Algorithms
+{
+/**
+ * Operator class performing Newton's iteration with standard step
+ * size control and adaptive matrix generation.
+ *
+ * This class performes a Newton iteration up to convergence
+ * determined by #control. If after an update the norm of the residual
+ * has become larger, then step size control is activated and the
+ * update is subsequently divided by two until the residual actually
+ * becomes smaller (or the minimal scaling factor determined by
+ * #n_stepsize_iterations is reached).
+ *
+ * Since assembling matrices, depending on the implementation, tends
+ * to be costly, this method applies an adaptive reassembling
+ * strategy. Only if the reduction factor for the residual is more
+ * than #threshold, the event #bad_derivative is submitted to
+ * #inverse_derivative. It is up to this object to implement
+ * reassembling accordingly.
+ *
+ * @author Guido Kanschat, 2006, 2010
+ */
+  template <class VECTOR>
+  class Newton : public Operator<VECTOR>
+  {
+    public:
+                                      /**
+                                       * Constructor, receiving the
+                                       * application computing the
+                                       * residual and solving the
+                                       * linear problem.
+                                       */
+      Newton (Operator<VECTOR>& residual, Operator<VECTOR>& inverse_derivative);
+      
+                                      /**
+                                       * Declare the parameter
+                                       * applicable to Newton's method.
+                                       */
+      void declare_parameters (ParameterHandler& param);
+      
+                                      /**
+                                       * Read the parameters.
+                                       */
+      void initialize (ParameterHandler& param);
+      
+                                      /**
+                                       * The actual Newton iteration.
+                                       */
+      virtual void operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in);
+      
+      virtual void notify(const Event&);
+
+                                      /**
+                                       * Set the maximal residual
+                                       * reduction allowed without
+                                       * triggering assembling in the
+                                       * next step. Return the
+                                       * previous value.
+                                       */
+      double threshold(double new_value);
+
+                                      /**
+                                       * Control object for the
+                                       * Newton iteration.
+                                       */
+      ReductionControl control;
+    private:
+                                      /**
+                                       * The operator computing the residual.
+                                       */
+      SmartPointer<Operator<VECTOR> > residual;
+
+                                      /**
+                                       * The operator applying the
+                                       * inverse derivative to the residual.
+                                       */
+      SmartPointer<Operator<VECTOR> > inverse_derivative;
+      
+                                      /**
+                                       * This flag is set by the
+                                       * function assemble(),
+                                       * indicating that the matrix
+                                       * must be assembled anew upon
+                                       * start.
+                                       */
+      bool assemble_now;
+
+                                      /**
+                                       * A flag used to decide how
+                                       * many stepsize iteration
+                                       * should be made. Default is
+                                       * the original value of 21.
+                                       *
+                                       * Enter zero here to turn of
+                                       * stepsize control.
+                                       *
+                                       * @note Controlled by
+                                       * <tt>Stepsize iterations</tt> in
+                                       * parameter file
+                                       */
+      unsigned int n_stepsize_iterations;
+      
+                                      /**
+                                       * Threshold for re-assembling matrix.
+                                       *
+                                       * If the quotient of two
+                                       * consecutive residuals is
+                                       * smaller than this threshold,
+                                       * the system matrix is not
+                                       * assembled in this step.
+                                       *
+                                       * @note This parameter should be
+                                       * adjusted to the residual gain
+                                       * of the inner solver.
+                                       *
+                                       * The default values is zero,
+                                       * resulting in reassembling in
+                                       * every Newton step.
+                                       */
+      double assemble_threshold;
+                                      /**
+                                       * Print residual, update and
+                                       * updated solution after each
+                                       * step into file
+                                       * <tt>Newton_NNN</tt>?
+                                       */
+      bool debug_vectors;
+                                      /**
+                                       * Write debug output to
+                                       * #deallog; the higher the
+                                       * number, the more output.
+                                       */
+      unsigned int debug;      
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/newton.templates.h b/deal.II/deal.II/include/numerics/newton.templates.h
new file mode 100644 (file)
index 0000000..1918473
--- /dev/null
@@ -0,0 +1,169 @@
+//---------------------------------------------------------------------------
+//    $Id: newton.h 753 2009-05-27 00:30:14Z kanschat $
+//
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/newton.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+  template <class VECTOR>
+  Newton<VECTOR>::Newton(Operator<VECTOR>& residual, Operator<VECTOR>& inverse_derivative)
+                 :
+                 residual(&residual), inverse_derivative(&inverse_derivative),
+                 assemble_now(false),
+                 n_stepsize_iterations(21),
+                 assemble_threshold(0.),
+                 debug_vectors(false),
+                 debug(0)
+  {}
+
+
+  template <class VECTOR>
+  void
+  Newton<VECTOR>::declare_parameters(ParameterHandler& param)
+  {
+    param.enter_subsection("Newton");
+    ReductionControl::declare_parameters (param);
+    param.declare_entry("Assemble threshold", "0.", Patterns::Double());
+    param.declare_entry("Stepsize iterations", "21", Patterns::Integer());
+    param.declare_entry("Debug level", "0", Patterns::Integer());
+    param.declare_entry("Debug vectors", "false", Patterns::Bool());
+    param.leave_subsection();    
+  }
+
+  template <class VECTOR>
+  void
+  Newton<VECTOR>::initialize (ParameterHandler& param)
+  {
+    param.enter_subsection("Newton");
+    control.parse_parameters (param);
+    assemble_threshold = param.get_double("Assemble threshold");
+    n_stepsize_iterations = param.get_double("Stepsize iterations");
+    debug_vectors = param.get_bool("Debug vectors");
+    param.leave_subsection ();
+  }
+  
+
+  template <class VECTOR>
+  void
+  Newton<VECTOR>::notify(const Event& e)
+  {
+    residual->notify(e);
+    inverse_derivative->notify(e);
+  }
+  
+  
+  template <class VECTOR>
+  void
+  Newton<VECTOR>::operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in)
+  {
+    Assert (out.size() == 1, ExcNotImplemented());
+    deallog.push ("Newton");
+
+    VECTOR& u = *out(0);
+  
+    if (debug>2)
+      deallog << "u: " << u.l2_norm() << std::endl;
+
+    GrowingVectorMemory<VECTOR> mem;
+    typename VectorMemory<VECTOR>::Pointer Du(mem);
+    typename VectorMemory<VECTOR>::Pointer res(mem);
+  
+    res->reinit(u);
+    NamedData<VECTOR*> src1;
+    NamedData<VECTOR*> src2;
+    VECTOR* p = &u;
+    src1.add(p, "Newton iterate");
+    src1.merge(in);
+    p = res;
+    src2.add(p, "Newton residual");
+    src2.merge(src1);
+    NamedData<VECTOR*> out1;
+    out1.add(p, "Residual");
+    p = Du;
+    NamedData<VECTOR*> out2;
+    out2.add(p, "Update");
+  
+    unsigned int step = 0;
+                                    // fill res with (f(u), v)
+    (*residual)(out1, src1);
+    double resnorm = res->l2_norm();
+    double old_residual = resnorm / assemble_threshold + 1;
+  
+    while (control.check(step++, resnorm) == SolverControl::iterate)
+      {
+                                        // assemble (Df(u), v)
+       if (resnorm/old_residual >= assemble_threshold)
+         inverse_derivative->notify (Events::bad_derivative);
+       
+       Du->reinit(u);
+       try 
+         {
+           (*inverse_derivative)(out2, src2);
+         }
+       catch (SolverControl::NoConvergence& e)
+         {
+           deallog << "Inner iteration failed after "
+                   << e.last_step << " steps with residual "
+                   << e.last_residual << std::endl;
+         }
+       
+       u.add(-1., *Du);
+       old_residual = resnorm;
+       (*residual)(out1, src1);
+       resnorm = res->l2_norm();
+       
+                                        // Step size control
+       unsigned int step_size = 0;
+       while (resnorm >= old_residual)
+         {
+           ++step_size;
+           if (step_size > n_stepsize_iterations)
+             {
+               deallog << "No smaller stepsize allowed!";
+               break;
+             }
+           if (control.log_history())
+             deallog << "Trying step size: 1/" << (1<<step_size)
+                     << " since residual was " << residual << std::endl;
+           u.add(1./(1<<step_size), *Du);
+           (*residual)(out1, src1);
+           resnorm = res->l2_norm();
+         }
+       
+       if (debug_vectors)
+         {
+           std::ostringstream streamOut;
+           streamOut << "Newton_" << std::setw(3) << std::setfill('0') << step;
+           std::string name = streamOut.str();
+           NamedData<VECTOR*> out;
+           out.add(&u, "solution");
+           out.add(Du, "update");
+           out.add(res, "residual");       
+//TODO: Implement!
+//         app->data_out(name, out);
+         }
+      }  
+    deallog.pop();
+
+                                    // in case of failure: throw
+                                    // exception
+    if (control.last_check() != SolverControl::success)
+      throw SolverControl::NoConvergence (control.last_step(),
+                                         control.last_value());
+                                    // otherwise exit as normal
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/include/numerics/operator.h b/deal.II/deal.II/include/numerics/operator.h
new file mode 100644 (file)
index 0000000..e6d70f7
--- /dev/null
@@ -0,0 +1,79 @@
+//---------------------------------------------------------------------------
+//    $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__operator_h
+#define __deal2__operator_h
+
+#include <base/config.h>
+#include <base/named_data.h>
+#include <numerics/event.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+/**
+ * The abstract base class of all algorithms in this library. An
+ * operator is an object with an #operator(), which transforms a set
+ * of named vectors into another set of named vectors.
+ *
+ * Furthermore, an operator can be notified of parameter changes by
+ * the calling routine. The outer iteration can notify() the Operator
+ * of an Event, which could be for instance a change of mesh, a
+ * different time step size or too slow convergence of Newton's
+ * method, which would then trigger reassembling of a matrix or
+ * similar things.
+ *
+ * @author Guido Kanschat, 2010
+ */
+  template <class VECTOR>
+  class Operator : public Subscriptor
+  {
+    public:
+                                      /**
+                                       * The virtual destructor.
+                                       */
+      ~Operator();
+
+                                      /**
+                                       * The actual operation, which
+                                       * is implemented in a derived class.
+                                       */
+      virtual void operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in) = 0;
+
+                                      /**
+                                       * Register an event triggered
+                                       * by an outer iteration.
+                                       */
+      virtual void notify(const Event&);
+                                      /**
+                                       * Clear all #notifications.
+                                       */
+      void clear_events();
+    protected:
+                                      /**
+                                       * Accumulate reasons for
+                                       * reassembling here. If any of
+                                       * those is set, the function
+                                       * solve() of a terminal
+                                       * application must take care
+                                       * of reassembling the matrix.
+                                       */
+      Event notifications;      
+  };
+
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/operator.templates.h b/deal.II/deal.II/include/numerics/operator.templates.h
new file mode 100644 (file)
index 0000000..fc9095b
--- /dev/null
@@ -0,0 +1,39 @@
+//---------------------------------------------------------------------------
+//    $Id: newton.h 753 2009-05-27 00:30:14Z kanschat $
+//
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/operator.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+  template <class VECTOR>
+  Operator<VECTOR>::~Operator()
+  {}
+  
+  
+  template <class VECTOR>
+  void Operator<VECTOR>::notify(const Event& e)
+  {
+    notifications += e;
+  }
+  
+  
+  template <class VECTOR>
+  void
+  Operator<VECTOR>::clear_events ()
+  {
+    notifications.clear();
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/event.all_dimensions.cc b/deal.II/deal.II/source/numerics/event.all_dimensions.cc
new file mode 100644 (file)
index 0000000..ad45b3a
--- /dev/null
@@ -0,0 +1,69 @@
+//---------------------------------------------------------------------------
+//    $Id: vectors.all_dimensions.cc 14038 2006-10-23 02:46:34Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/event.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+//TODO: Thread safety
+
+namespace Algorithms
+{
+  std::vector<std::string> Event::names;
+
+  Event
+  Event::assign(const char* name)
+  {
+    unsigned int index = names.size();
+    names.push_back(name);
+
+    Event result;
+                                    // The constructor generated an
+                                    // object with all flags equal
+                                    // zero. Now we set the new one.
+    result.flags[index] = true;
+
+    return result;
+  }
+
+
+  Event::Event ()
+                 :
+                 all_true(false),
+                 flags(names.size(), false)
+  {}
+
+
+  void
+  Event::clear ()
+  {
+    all_true = false;
+    std::fill(flags.begin(), flags.end(), false);
+  }
+
+
+  void
+  Event::all ()
+  {
+    all_true = true;
+  }
+
+  namespace Events
+  {
+    const Event bad_derivative = Event::assign("Bad Derivative");
+    const Event new_time = Event::assign("New Time");
+    const Event new_time_step_size = Event::assign("New Time Step Size");
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
index a89e80acc964cc3952b86c44117ee90f3f424218..1ccad5b062677fd13d413bb2fe5cc34cca982bc7 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
diff --git a/deal.II/deal.II/source/numerics/operators.cc b/deal.II/deal.II/source/numerics/operators.cc
new file mode 100644 (file)
index 0000000..0e09fb4
--- /dev/null
@@ -0,0 +1,30 @@
+//---------------------------------------------------------------------------
+//    $Id: vectors.cc 18474 2009-03-10 16:55:43Z kronbichler $
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <base/parameter_handler.h>
+#include <base/logstream.h>
+#include <lac/vector_memory.h>
+
+#include <numerics/operator.templates.h>
+#include <numerics/newton.templates.h>
+
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+using namespace Algorithms;
+
+#include "operators.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/operators.inst.in b/deal.II/deal.II/source/numerics/operators.inst.in
new file mode 100644 (file)
index 0000000..75979cd
--- /dev/null
@@ -0,0 +1,18 @@
+//---------------------------------------------------------------------------
+//    $Id: vectors.inst.in 19046 2009-07-08 19:30:23Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+for (VEC : SERIAL_VECTORS)
+{  
+  template class Operator<VEC>;
+  template class Newton<VEC>;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.