]> https://gitweb.dealii.org/ - dealii.git/commitdiff
working on new example
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 5 Jul 2010 18:14:15 +0000 (18:14 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 5 Jul 2010 18:14:15 +0000 (18:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@21454 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/doxygen/Makefile
deal.II/examples/doxygen/theta_timestepping.cc
deal.II/examples/step-20/doc/intro.dox
deal.II/examples/step-20/step-20.cc

index 1ed195f86c1fc543d1898c1fd819812aa1e14567..5bd22b1a7ca00b49be1bcfbab12986f0ff6bea2f 100644 (file)
@@ -2,7 +2,7 @@
 #    $Id$
 #    Version: $Name$
 #
-#    Copyright (C) 1998 - 2007 by the deal.II authors
+#    Copyright (C) 1998 - 2007, 2010 by the deal.II authors
 #
 #    This file is subject to QPL and may not be  distributed
 #    without copyright and license information. Please refer
@@ -13,7 +13,8 @@
 include ../../common/Make.global_options
 
 all: product_matrix$(EXEEXT) block_matrix_array$(EXEEXT) \
-     compressed_block_sparsity_pattern$(EXEEXT)
+     compressed_block_sparsity_pattern$(EXEEXT) \
+     theta_timestepping$(EXEEXT)
 
 ######################################################################
 # Compilation of source code
@@ -42,6 +43,10 @@ compressed_block_sparsity_pattern$(EXEEXT): compressed_block_sparsity_pattern.g.
        @echo ============================ Linking $@
        @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
 
+theta_timestepping$(EXEEXT): theta_timestepping.g.$(OBJEXT) $(lib-deal2-2d.g) $(lib-lac.g) $(lib-base.g)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
 ######################################################################
 # Pseudo target for cleaning directory
 ######################################################################
index 8d9bf626aa36b430dc106bee9ae4e9beae422ada..6f42cbd37fac5f872383c1d66deee52657f19c49 100644 (file)
 using namespace dealii;
 using namespace Algorithms;
 
+
 class Explicit
-  : Operator<Vector<double> >
+  : 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&);
+    void initialize_timestep_data(const TimestepData&);
   private:
-    const TimeStepData* timestep_data;
+    const TimestepData* timestep_data;
+    SmartPointer<const FullMatrix<double>, Explicit> matrix;
+    FullMatrix<double> m;
 };
 
   
 class Implicit
-  : Operator<Vector<double> >
+  : 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&);
+    void initialize_timestep_data(const TimestepData&);
   private:
-    const TimeStepData* timestep_data;
+    const TimestepData* timestep_data;
+    SmartPointer<const FullMatrix<double>, Implicit> matrix;
+    FullMatrix<double> m;
 };
 
   
 int main()
 {
-  Explicit op_explicit;
-  Implicit op_implicit;
+  FullMatrix<double> matrix(2);
+  matrix(0,0) = -.1;
+  matrix(1,1) = -.1;
+  matrix(0,1) = 31.4;
+  matrix(1,0) = -31.4;
+
+  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.notify(Events::initial);
+  
   Vector<double> value(2);
   NamedData<Vector<double>*> indata;
   NamedData<Vector<double>*> outdata;
-  outdata.add(&value, "value");
+  Vector<double>* p = &value;
+  outdata.add(p, "value");
   solver(outdata, indata);
 }
+
+
+Explicit::Explicit(const FullMatrix<double>& M)
+               :
+               matrix(&M)
+{}
+
+
+void
+Explicit::initialize_timestep_data(const TimestepData& t)
+{
+  timestep_data = &t;
+}
+
+
+void
+Explicit::operator() (NamedData<Vector<double>*>& out, const NamedData<Vector<double>*>& in)
+{
+  this->notifications.print(deallog);
+  deallog << std::endl;
+  unsigned int i = in.find("Previous time");
+  m.vmult(*out(0), *in(i));
+}
+
+
+Implicit::Implicit(const FullMatrix<double>& M)
+               :
+               matrix(&M)
+{}
+
+
+void
+Implicit::initialize_timestep_data(const TimestepData& t)
+{
+  timestep_data = &t;
+}
+
+
+void
+Implicit::operator() (NamedData<Vector<double>*>& out, const NamedData<Vector<double>*>& in)
+{
+  this->notifications.print(deallog);
+  deallog << std::endl;
+  unsigned int i = in.find("Previous time");
+  m.vmult(*out(0), *in(i));
+}
+
+
index 91c7b3ca77ae3127560d9f6ee4b340f16d30350d..28e6b48b88160580e029edc4cacd9726efa5b037 100644 (file)
@@ -525,43 +525,6 @@ this looks almost as expensive as solving with $S$ right away. However, note
 that in the inner iteration, we do not have to calculate $M^{-1}$, but only
 the inverse of its diagonal, which is cheap.
 
-To implement something like this, let us first generalize the
-<code>InverseMatrix</code> class so that it can work not only with
-<code>SparseMatrix</code> objects, but with any matrix type. This looks like so:
-
-@code
-template <class Matrix>
-class InverseMatrix
-{
-  public:
-    InverseMatrix (const Matrix &m);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const Matrix> matrix;
-
-    //...
-};
-
-
-template <class Matrix>
-void InverseMatrix<Matrix>::vmult (Vector<double>       &dst,
-                                   const Vector<double> &src) const
-{
-  SolverControl solver_control (src.size(), 1e-8*src.l2_norm());
-  SolverCG<> cg (solver_control, vector_memory);
-
-  dst = 0;
-  
-  cg.solve (*matrix, dst, src, PreconditionIdentity());        
-}
-@endcode
-
-Essentially, the only change we have made is the introduction of a template
-argument that generalizes the use of <code>SparseMatrix</code>.
-
 The next step is to define a class that represents the approximate Schur
 complement. This should look very much like the Schur complement class itself,
 except that it doesn't need the object representing $M^{-1}$ any more:
@@ -662,9 +625,7 @@ chose to introduce this material here anyway to demonstrate how to work with
 block matrices and to develop solvers and preconditioners, rather than using
 black box components from the library.
 
-For those interested in looking up the corresponding library classes: the
-<code>InverseMatrix</code> is roughly equivalent to the
-<code>PreconditionLACSolver</code> class in the library. Likewise, the Schur
+For those interested in looking up the corresponding library classes: the Schur
 complement class corresponds to the <code>SchurMatrix</code> class.
 
 
index febeb0602f06d83387a3efc2e1cf8bf9f0df948b..6cc260a81bddc9307c8abe8ecf8cdff5cb9f6a8a 100644 (file)
@@ -34,6 +34,7 @@
                                 // inverse of a matrix by calling an
                                 // iterative solver.
 #include <lac/iterative_inverse.h>
+#include <lac/schur_matrix.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
@@ -726,86 +727,6 @@ void MixedLaplaceProblem<dim>::assemble_system ()
                                  // rather only comment on
                                  // implementational aspects.
 
-
-                                 // @sect4{The <code>SchurComplement</code> class template}
-
-                                 // The next class is the Schur
-                                 // complement class. Its rationale
-                                 // has also been discussed in length
-                                 // in the introduction. The only
-                                 // things we would like to note is
-                                 // that the class, too, is derived
-                                 // from the <code>Subscriptor</code> class and
-                                 // that as mentioned above it stores
-                                 // pointers to the entire block
-                                 // matrix and the inverse of the mass
-                                 // matrix block using
-                                 // <code>SmartPointer</code> objects.
-                                 //
-                                 // The <code>vmult</code> function requires
-                                 // two temporary vectors that we do
-                                 // not want to re-allocate and free
-                                 // every time we call this
-                                 // function. Since here, we have full
-                                 // control over the use of these
-                                 // vectors (unlike above, where a
-                                 // class called by the <code>vmult</code>
-                                 // function required these vectors,
-                                 // not the <code>vmult</code> function
-                                 // itself), we allocate them
-                                 // directly, rather than going
-                                 // through the <code>VectorMemory</code>
-                                 // mechanism. However, again, these
-                                 // member variables do not carry any
-                                 // state between successive calls to
-                                 // the member functions of this class
-                                 // (i.e., we never care what values
-                                 // they were set to the last time a
-                                 // member function was called), we
-                                 // mark these vectors as <code>mutable</code>.
-                                 //
-                                 // The rest of the (short)
-                                 // implementation of this class is
-                                 // straightforward if you know the
-                                 // order of matrix-vector
-                                 // multiplications performed by the
-                                 // <code>vmult</code> function:
-class SchurComplement : public Subscriptor
-{
-  public:
-    SchurComplement (const BlockSparseMatrix<double> &A,
-                     const IterativeInverse<Vector<double> > &Minv);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
-    
-    mutable Vector<double> tmp1, tmp2;
-};
-
-
-SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
-                                  const IterativeInverse<Vector<double> > &Minv)
-                :
-                system_matrix (&A),
-                m_inverse (&Minv),
-                tmp1 (A.block(0,0).m()),
-                tmp2 (A.block(0,0).m())
-{}
-
-
-void SchurComplement::vmult (Vector<double>       &dst,
-                             const Vector<double> &src) const
-{
-  system_matrix->block(0,1).vmult (tmp1, src);
-  m_inverse->vmult (tmp2, tmp1);
-  system_matrix->block(1,0).vmult (dst, tmp2);
-}
-
-
                                  // @sect4{The <code>ApproximateSchurComplement</code> class template}
 
                                  // The third component of our solver
@@ -901,7 +822,11 @@ void MixedLaplaceProblem<dim>::solve ()
   m_inverse.solver.select("cg");
   ReductionControl inner_control(1000, 0., 1.e-13);
   m_inverse.solver.control = inner_control;
-  
+
+  SchurComplement schur_complement (m_inverse,
+                                   system_matrix.block(1,0),
+                                   system_matrix.block(1,0));
+
   Vector<double> tmp (solution.block(0).size());
 
                                    // Now on to the first
@@ -925,8 +850,6 @@ void MixedLaplaceProblem<dim>::solve ()
     schur_rhs -= system_rhs.block(1);
 
     
-    SchurComplement
-      schur_complement (system_matrix, m_inverse);
     
     ApproximateSchurComplement
       approximate_schur_complement (system_matrix);

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.