]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
more multigrid
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 30 Apr 1999 14:09:21 +0000 (14:09 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 30 Apr 1999 14:09:21 +0000 (14:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@1236 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/doc/parser/dealparse.y
deal.II/lac/include/lac/mgbase.h
deal.II/lac/source/mgbase.cc
tests/deal.II/mg2.cc [new file with mode: 0644]
tests/lac/Makefile
tests/lac/mg.cc [new file with mode: 0644]
tests/lac/mgbase.cc [moved from tests/deal.II/mg1.cc with 66% similarity]
tests/lac/solver.cc
tests/lac/testmatrix.cc
tests/lac/testmatrix.h

index e864500b26fc6c0e8a9950cff2ea8591c356bbc1..e0e83596d5f06da041c91d03d6d6f17b4cec62f8 100644 (file)
@@ -131,7 +131,7 @@ template<class MATRIX>
 class MGMatrix
   :
   public Subscriptor,
-  public Vector<MATRIX>
+  public vector<MATRIX>
 {
   public:
                                     /**
@@ -159,7 +159,7 @@ class MGCoarseGridLACIteration
                                      * preconditioning method for later
                                      * use.
                                      */
-    MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&);
+    MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&);
     
                                     /**
                                      * Implementation of the abstract
@@ -174,7 +174,7 @@ class MGCoarseGridLACIteration
                                     /**
                                      * Reference to the solver.
                                      */
-    const SOLVER& solver;
+    SOLVER& solver;
                                     /**
                                      * Reference to the matrix.
                                      */
@@ -280,7 +280,7 @@ class MGBase
 
 template<class SOLVER, class MATRIX, class PRECOND>
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
-::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p)
+::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
                :
                solver(s), matrix(m), precondition(p)
 {}
@@ -288,7 +288,7 @@ MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
 template<class SOLVER, class MATRIX, class PRECOND>
 void
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int level, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<float>& dst, const Vector<float>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
index 1ccc5258a6dd2f09d11e5436e87b5b3f7fdd9bf6..c6252fa0d8d1dafe6c3f5ecff1b52a573581a262 100644 (file)
@@ -45,3 +45,8 @@ MGBase::level_mgstep(unsigned level,
   post_smooth.smooth(level, s[level], d[level]);
 }
 
+MGTransferBase::~MGTransferBase()
+{}
+
+MGSmootherBase::~MGSmootherBase()
+{}
index c8d0bafd2dcf9f3662dd759a2b6107f9dfc1ddab..e84390d1d5d141cb0a86c236d2c62ad9723d2df8 100644 (file)
@@ -197,6 +197,8 @@ enumerator: IDENTIFIER
 class_declaration: class_head
   | class_head
     {
+      cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl
+          << setw(OUTW) << "@In-Class:" << class_stack.top() << endl;
       class_name = class_stack.top() + string("::") + class_name;
       class_stack.push(class_name);
       access_stack.push(string("private"));
@@ -204,8 +206,6 @@ class_declaration: class_head
     {
       class_stack.pop();
       access_stack.pop();
-      cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl
-          << setw(OUTW) << "@In-Class:" << class_stack.top() << endl;
     }
 ;
 
index f614a308dcc88bae978da27ef1aee1053fad7880..3cb773a6b2a738a97771f6dc8304a26bb013aa97 100644 (file)
@@ -131,7 +131,7 @@ template<class MATRIX>
 class MGMatrix
   :
   public Subscriptor,
-  public Vector<MATRIX>
+  public vector<MATRIX>
 {
   public:
                                     /**
@@ -159,7 +159,7 @@ class MGCoarseGridLACIteration
                                      * preconditioning method for later
                                      * use.
                                      */
-    MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&);
+    MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&);
     
                                     /**
                                      * Implementation of the abstract
@@ -174,7 +174,7 @@ class MGCoarseGridLACIteration
                                     /**
                                      * Reference to the solver.
                                      */
-    const SOLVER& solver;
+    SOLVER& solver;
                                     /**
                                      * Reference to the matrix.
                                      */
@@ -280,7 +280,7 @@ class MGBase
 
 template<class SOLVER, class MATRIX, class PRECOND>
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
-::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p)
+::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
                :
                solver(s), matrix(m), precondition(p)
 {}
@@ -288,7 +288,7 @@ MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
 template<class SOLVER, class MATRIX, class PRECOND>
 void
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int level, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<float>& dst, const Vector<float>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
index 1ccc5258a6dd2f09d11e5436e87b5b3f7fdd9bf6..c6252fa0d8d1dafe6c3f5ecff1b52a573581a262 100644 (file)
@@ -45,3 +45,8 @@ MGBase::level_mgstep(unsigned level,
   post_smooth.smooth(level, s[level], d[level]);
 }
 
+MGTransferBase::~MGTransferBase()
+{}
+
+MGSmootherBase::~MGSmootherBase()
+{}
diff --git a/tests/deal.II/mg2.cc b/tests/deal.II/mg2.cc
new file mode 100644 (file)
index 0000000..6e3ace4
--- /dev/null
@@ -0,0 +1,22 @@
+// $Id$
+
+// deal_II_libraries.g=-ldeal_II_2d.g
+// deal_II_libraries=-ldeal_II_2d
+
+#include <numerics/multigrid.h>
+#include <numerics/mg_smoother.h>
+#include <grid/tria.h>
+#include <grid/mg_dof.h>
+#include <fe/fe_lib.lagrange.h>
+
+main()
+{
+  Triangulation<2> tr;
+  MGDoFHandler<2> dof(&tr);
+  FELinear<2> fe;
+  
+  tr.create_hypercube(-1.,1.);
+  tr.refine_global(3);
+  dof.distribute_dofs(fe);
+}
+
index 3d1a1cf7abe0de9de533ea06e0c9173631bbb230..1b7ce6c299adb20814f967a6192d40ee25c09c34 100644 (file)
@@ -41,7 +41,7 @@ libraries = $(libs)
 flags     = $(CXXFLAGS.o) $(CXXFLAGS)
 endif
 
-all: solver
+all: solver mgbase mg
 
 ############################################################
 # Typical block for building a running program
@@ -67,6 +67,39 @@ endif
 solver: $(solver-o-files) $(libraries)
        $(CXX) $(flags) -o $@ $^
 
+
+############################################################
+
+
+mgbase-cc-files = mgbase.cc
+
+ifeq ($(debug-mode),on)
+mgbase-o-files = $(mgbase-cc-files:.cc=.go)
+else
+mgbase-o-files = $(mgbase-cc-files:.cc=.o)
+endif
+
+mgbase: $(mgbase-o-files) $(libraries)
+       $(CXX) $(flags) -o $@ $^
+
+
+
+############################################################
+
+
+mg-cc-files = mg.cc testmatrix.cc
+
+ifeq ($(debug-mode),on)
+mg-o-files = $(mg-cc-files:.cc=.go)
+else
+mg-o-files = $(mg-cc-files:.cc=.o)
+endif
+
+mg: $(mg-o-files) $(libraries)
+       $(CXX) $(flags) -o $@ $^
+
+
+
 ############################################################
 # Continue with other targets if needed
 ############################################################
@@ -99,7 +132,7 @@ veryclean: clean
 # Automatic generation of dependencies
 ############################################################
 
-all-cc-files = $(solver-cc-files) # $(target2-cc-files) ...
+all-cc-files = $(solver-cc-files) $(mgbase-cc-files) $(mg-cc-files) # $(target2-cc-files) ...
 
 Make.depend: $(all-cc-files)
        @echo =====Dependecies=== Make.depend
diff --git a/tests/lac/mg.cc b/tests/lac/mg.cc
new file mode 100644 (file)
index 0000000..439d47f
--- /dev/null
@@ -0,0 +1,160 @@
+// $Id$
+
+// Test program for multigrid with test matrices (fd)
+
+#include <cmath>
+#include <fstream>
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparsematrix.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+#include <lac/solver_control.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/solver_richardson.h>
+#include <lac/precondition.h>
+#include <lac/mgbase.h>
+
+class FDMG
+  :
+  public MGBase
+{
+  private:
+    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+  public:
+    FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+        FDMGTransfer& transfer)
+                   :
+                   MGBase(transfer, maxlevel, 0),
+                   matrices(&matrices)
+      {}
+    
+    
+    virtual void level_residual (unsigned int level,
+                                Vector<float>& dst,
+                                const Vector<float>& src,
+                                const Vector<float>& rhs);
+};
+
+class MGSmootherLAC
+  :
+  MGSmootherBase
+{
+  private:
+    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+  public:
+    MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+    
+    virtual void smooth (const unsigned int level,
+                        Vector<float> &u,
+                        const Vector<float> &rhs) const;
+    
+};
+
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float>  >,
+SparseMatrix<float>, PreconditionRelaxation<SparseMatrix<float> , Vector<float> > >
+Coarse;
+
+class MGPrecondition
+{
+  private:
+    FDMG& mg;
+    MGSmootherLAC& smooth;
+    Coarse& coarse;
+  public:
+    MGPrecondition(FDMG& m, MGSmootherLAC& s, Coarse& c)
+                   :
+                   mg(m), smooth(s), coarse(c)
+      {}
+    
+    
+    void operator () (Vector<double>& dst, const Vector<double>& src) const
+      {
+       
+      }
+    
+    
+};
+
+
+main()
+{
+  ofstream logfile("mg.output");
+  deallog.attach(logfile);
+  
+  PrimitiveVectorMemory<Vector<double>  > mem;
+  SolverControl control(100, 1.e-5, true);
+
+  const unsigned int base = 3;
+  const unsigned int maxlevel = 6;
+  
+  const unsigned int maxsize = base * (1 << maxlevel);
+  
+                                  // grid transfer
+  FDMGTransfer transfer(maxsize, maxsize, maxlevel);
+
+                                  // coarse grid solver
+  PrimitiveVectorMemory<Vector<float> > cgmem;
+  SolverControl cgcontrol(100, 1.e-12);
+  SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+
+  vector<SparseMatrixStruct >  structure(maxlevel+1);
+  MGMatrix<SparseMatrix<float> > A(maxlevel+1);
+  
+  for (unsigned int level = 0; level <= maxlevel; ++level)
+    {
+      const unsigned int size = base * (1 << level);
+
+      const unsigned int dim = (size+1)*(size+1);
+
+                                      // Make matrix
+      FDMatrix testproblem(size, size);
+      structure[level].reinit(dim, dim, 5);
+      testproblem.build_structure(structure[level]);
+      structure[level].compress();
+      A[level].reinit(structure[level]);
+      testproblem.laplacian(A[level]);
+
+      FDMG multigrid(level, A, transfer);
+
+      PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
+       cgprec(A[0], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+      
+      Coarse coarsegrid(cgcg, A[0], cgprec);
+      
+      MGSmootherLAC smoother(A);
+      
+//      MGPrecondition prec(multigrid);
+    }
+}
+
+void
+FDMG::level_residual (unsigned int level,
+                     Vector<float>& dst,
+                     const Vector<float>& src,
+                     const Vector<float>& rhs)
+{
+  (*matrices)[level].residual(dst, src, rhs);
+}
+
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+               :
+               matrices(&matrix)
+{}
+
+
+void
+MGSmootherLAC::smooth (const unsigned int level,
+                      Vector<float> &u,
+                      const Vector<float> &rhs) const
+{
+  SolverControl control(0.,1);
+  PrimitiveVectorMemory<Vector<float> > mem;
+  SolverRichardson<SparseMatrix<float> , Vector<float>  > rich(control, mem);
+  PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
+    prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+
+  rich.solve((*matrices)[level], u, rhs, prec);
+}
similarity index 66%
rename from tests/deal.II/mg1.cc
rename to tests/lac/mgbase.cc
index 8c2efd4f32c42fbeb5085ef3291a7bae75969e7b..acd703a6c856047893e31e44c12211220f92dc47 100644 (file)
@@ -3,8 +3,7 @@
 // deal_II_libraries.g=-ldeal_II_2d.g
 // deal_II_libraries=-ldeal_II_2d
 
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
+#include <lac/mgbase.h>
 
 class TransferTest
   :
@@ -16,7 +15,7 @@ class TransferTest
     virtual void restrict(unsigned l,
                          Vector<float>& dst,
                          const Vector<float>& src) const;
-    friend class MultiGridBase;
+    friend class MGBase;
 };
 
 class SmoothTest
@@ -29,17 +28,29 @@ class SmoothTest
                             const Vector<float>& rhs) const;
 };
 
+class MGCoarseGridTest
+  :
+  public MGCoarseGridSolver
+{
+  public:
+    virtual void operator () (unsigned int level,
+                             Vector<float>& dst,
+                             const Vector<float>& src) const ;
+};
+
 
 class MGTest
   :
-  public MultiGridBase
+  public MGBase
 {
     MGSmootherBase& smoother;
-    
+    MGCoarseGridSolver& solver;
+
   public:
-    MGTest(TransferTest& t, SmoothTest& s)
-                   : MultiGridBase(t, 9, 3),
-                     smoother(s)
+    MGTest(TransferTest& t, SmoothTest& s, MGCoarseGridTest& c)
+                   : MGBase(t, 9, 3),
+                     smoother(s),
+                     solver(c)
       {}
     
     virtual void level_residual(unsigned level,
@@ -49,7 +60,7 @@ class MGTest
     
     void doit()
       {
-       level_mgstep(9, smoother, smoother);
+       level_mgstep(9, smoother, smoother, solver);
       }
 };
 
@@ -57,7 +68,9 @@ main()
 {
   TransferTest tr;
   SmoothTest s;
-  MGTest mg(tr, s);
+  MGCoarseGridTest c;
+  
+  MGTest mg(tr, s, c);
   mg.doit();
   
 }
@@ -83,7 +96,7 @@ SmoothTest::smooth (const unsigned int  level,
                    Vector<float>      &,
                    const Vector<float>&) const
 {
-  cout << "Smoothing on" << level << endl;
+  cout << "Smoothing on " << level << endl;
 }
 
 void
@@ -94,3 +107,11 @@ MGTest::level_residual(unsigned l,
 {
   cout << "Residual on  " << l << endl;
 }
+
+void
+MGCoarseGridTest::operator() (unsigned int level,
+                             Vector<float>&,
+                             const Vector<float>&) const
+{
+  cout << "Solving on   " << level << endl;
+}
index a70a6c478749617e26ba805c2fb3515d9f5ec567..acaa7ed9b83accea4285fe16e31780f08d0bc36c 100644 (file)
@@ -3,6 +3,7 @@
 // Test program for linear solvers.
 
 #include <cmath>
+#include <fstream>
 #include "testmatrix.h"
 #include <base/logstream.h>
 #include <lac/sparsematrix.h>
 #include <lac/solver_richardson.h>
 #include <lac/precondition.h>
 
-#define MATRIX SparseMatrix<float> 
-#define VECTOR Vector<double> 
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_method( SOLVER& solver, const MATRIX& A,
+            VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+  u = 0.;
+  f = 1.;
+  solver.solve(A,u,f,P);
+}
 
 main()
 {
-  PrimitiveVectorMemory<VECTOR > mem;
+  ofstream logfile("solver.output");
+  deallog.attach(logfile);
+  
+  PrimitiveVectorMemory<Vector<double>  > mem;
   SolverControl control(100, 1.e-5, true);
-  SolverCG<MATRIX, VECTOR > cg(control, mem);
-  SolverGMRES<MATRIX, VECTOR > gmres(control, mem,20);
-  SolverBicgstab<MATRIX, VECTOR > bicgstab(control, mem);
-  SolverRichardson<MATRIX, VECTOR > rich(control, mem);
+  SolverCG<SparseMatrix<float> , Vector<double>  > cg(control, mem);
+  SolverGMRES<SparseMatrix<float> , Vector<double>  > gmres(control, mem,20);
+  SolverBicgstab<SparseMatrix<float> , Vector<double>  > bicgstab(control, mem);
+  SolverRichardson<SparseMatrix<float> , Vector<double>  > rich(control, mem);
 
   for (unsigned int size=10; size <= 10; size *= 10)
     {
@@ -38,49 +49,40 @@ main()
       SparseMatrixStruct structure(dim, dim, 5);
       testproblem.build_structure(structure);
       structure.compress();
-      MATRIX A(structure);
+      SparseMatrix<float>  A(structure);
       testproblem.laplacian(A);
 
-      PreconditionIdentity<VECTOR >
+      PreconditionIdentity<Vector<double>  >
        prec_no;
-      PreconditionRelaxation<MATRIX, VECTOR>
-       prec_ssor(A, &MATRIX::template precondition_SSOR<double>, 1.2);
+      PreconditionRelaxation<SparseMatrix<float> , Vector<double> >
+       prec_sor(A, &SparseMatrix<float> ::template precondition_SOR<double>, 1.2);
+      PreconditionRelaxation<SparseMatrix<float> , Vector<double> >
+       prec_ssor(A, &SparseMatrix<float> ::template precondition_SSOR<double>, 1.2);
       
-      VECTOR f(dim);
-      VECTOR u(dim);
+      Vector<double>  f(dim);
+      Vector<double>  u(dim);
       
       deallog.push("no");
 
-      f = 1.;
-      u = 0.;
-      cg.solve(A,u,f,prec_no);
-
-      f = 1.;
-      u = 0.;
-      bicgstab.solve(A,u,f,prec_no);
-
-      f = 1.;
-      u = 0.;
-      gmres.solve(A,u,f,prec_no);
+      check_method(cg,A,u,f,prec_no);
+      check_method(bicgstab,A,u,f,prec_no);
+      check_method(gmres,A,u,f,prec_no);
 
       deallog.pop();
       deallog.push("ssor");      
 
-      f = 1.;
-      u = 0.;
-      rich.solve(A,u,f,prec_ssor);
-
-      f = 1.;
-      u = 0.;
-      cg.solve(A,u,f,prec_ssor);
+      check_method(rich,A,u,f,prec_ssor);
+      check_method(cg,A,u,f,prec_ssor);
+      check_method(bicgstab,A,u,f,prec_ssor);
+      check_method(gmres,A,u,f,prec_ssor);
 
-      f = 1.;
-      u = 0.;
-      bicgstab.solve(A,u,f,prec_ssor);
+      deallog.pop();
+      deallog.push("sor");      
 
-      f = 1.;
-      u = 0.;
-      gmres.solve(A,u,f,prec_ssor);
+      check_method(rich,A,u,f,prec_sor);
+      check_method(cg,A,u,f,prec_sor);
+      check_method(bicgstab,A,u,f,prec_sor);
+      check_method(gmres,A,u,f,prec_sor);
 
       deallog.pop();
     }
index 656cadef33095331dc444bda91cc25a774a067c5..8b22d664e731ffc21bde3c2c1117d18d490e6120 100644 (file)
@@ -78,5 +78,142 @@ FDMatrix::laplacian(SparseMatrix<number>& A) const
     } 
 }
 
+FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
+                          unsigned int nlevels)
+               :
+               structures(nlevels-1), matrices(nlevels-1)
+{
+  unsigned int power = 1 << nlevels;
+  
+  Assert ((nx%power)==0, ExcDivide(nx,power));
+  Assert ((ny%power)==0, ExcDivide(ny,power));
+  
+  nx /= power;
+  ny /= power;
+  
+  for (unsigned int level = 0; level < nlevels-1; ++ level)
+    {
+      build_matrix(nx,ny,structures[level],matrices[level]);
+      nx *= 2;
+      ny *= 2;
+    }
+}
+
+void
+FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
+                          SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+{
+  structure.reinit((nx+1)*(ny+1),(2*nx+1)*(2*ny+1),9);
+  
+                                  // Traverse all points of coarse grid
+  for (unsigned int i=0 ; i<=ny; ++i)
+    for (unsigned int j=0 ; j<=nx; ++j)
+      {
+                                        // coarse grid point number
+       unsigned int ncoarse =j+(nx+1)*i;
+                                        // same point on fine grid
+       unsigned int nfine = 2*j+(4*nx+2)*i;
+       
+       structure.add(ncoarse,nfine);
+                                        // left
+       if (j>0)
+         {
+           structure.add(ncoarse,nfine-1);
+                                            // lower left
+           if (i>0)
+             structure.add(ncoarse,nfine-(2*nx+1)-1);
+                                            // upper left
+           if (i<ny)
+             structure.add(ncoarse,nfine+(2*nx+1)-1);
+         }
+                                        // right
+       if (j<nx)
+         {
+           structure.add(ncoarse, nfine+1);
+                                            // lower right
+           if (i>0)
+             structure.add(ncoarse,nfine-(2*nx+1)+1);
+                                            // upper right
+           if (i<ny)
+             structure.add(ncoarse,nfine+(2*nx+1)+1);
+         }
+
+                                            // lower
+           if (i>0)
+             structure.add(ncoarse,nfine-(2*nx+1));
+                                            // upper
+           if (i<ny)
+             structure.add(ncoarse,nfine+(2*nx+1));
+      }
+
+  structure.compress();
+  matrix.reinit(structure);
+  
+  for (unsigned int i=0 ; i<=ny; ++i)
+    for (unsigned int j=0 ; j<=nx; ++j)
+      {
+                                        // coarse grid point number
+       unsigned int ncoarse =j+(nx+1)*i;
+                                        // same point on fine grid
+       unsigned int nfine = 2*j+(4*nx+2)*i;
+       
+       matrix.set(ncoarse,nfine,1.);
+                                        // left
+       if (j>0)
+         {
+           matrix.set(ncoarse,nfine-1,.5);
+                                            // lower left
+           if (i>0)
+             matrix.set(ncoarse,nfine-(2*nx+1)-1,.25);
+                                            // upper left
+           if (i<ny)
+             matrix.set(ncoarse,nfine+(2*nx+1)-1,.25);
+         }
+                                        // right
+       if (j<nx)
+         {
+           matrix.set(ncoarse,nfine+1,.5);
+                                            // lower right
+           if (i>0)
+             matrix.set(ncoarse,nfine-(2*nx+1)+1,.25);
+                                            // upper right
+           if (i<ny)
+             matrix.set(ncoarse,nfine+(2*nx+1)+1,.25);
+         }
+
+                                            // lower
+           if (i>0)
+             matrix.set(ncoarse,nfine-(2*nx+1),.5);
+                                            // upper
+           if (i<ny)
+             matrix.set(ncoarse,nfine+(2*nx+1),.5);
+      }
+  
+}
+
+void
+FDMGTransfer::prolongate (const unsigned int   to_level,
+                         Vector<float>       &dst,
+                         const Vector<float> &src) const
+{
+  Assert((to_level>0) && (to_level<=matrices.size()),
+        ExcIndexRange(to_level, 0, matrices.size()+1));
+  
+  matrices[to_level-1].Tvmult(dst,src);
+}
+
+
+void
+FDMGTransfer::restrict (const unsigned int   from_level,
+                       Vector<float>       &dst,
+                       const Vector<float> &src) const
+{
+  Assert((from_level>0) && (from_level<=matrices.size()),
+        ExcIndexRange(from_level, 0, matrices.size()+1));
+
+  matrices[from_level-1].vmult(dst,src);
+}
+
+
 template void FDMatrix::laplacian(SparseMatrix<double>&) const;
 template void FDMatrix::laplacian(SparseMatrix<float>&) const;
index ffb32d8e8a3bbc89d806d5e9727fb6ee2d374cf9..641f949e0fa0438d04f9b50a5b84d703dad99a70 100644 (file)
@@ -1,6 +1,8 @@
 // $Id$
 
 #include <lac/forward-declarations.h>
+#include <base/exceptions.h>
+#include <lac/mgbase.h>
 
 /**
  * Finite difference matrix on uniform grid.
@@ -37,3 +39,64 @@ class FDMatrix
                                      */
     unsigned int ny;
 };
+
+/**
+ * Grid transfer for finite differences on uniform grid.
+ */
+
+class FDMGTransfer
+  :
+  public MGTransferBase
+{
+  public:
+                                    /**
+                                     * Constructor. Prepares grid
+                                     * transfer matrices for #nlevels#
+                                     * levels on an #nx# times #ny#
+                                     * grid.
+                                     */
+    FDMGTransfer(unsigned int nx, unsigned int ny,
+                unsigned int nlevels);
+
+                                    /**
+                                     * Implementation of abstract
+                                     * function in #MGTranferBase#.
+                                     */
+    virtual void prolongate (const unsigned int   to_level,
+                            Vector<float>       &dst,
+                            const Vector<float> &src) const;
+
+                                    /**
+                                     * Implementation of abstract
+                                     * function in #MGTranferBase#.
+                                     */
+    virtual void restrict (const unsigned int   from_level,
+                          Vector<float>       &dst,
+                          const Vector<float> &src) const;
+
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException2(ExcDivide, unsigned int, unsigned int,
+                  << "Cannot divide " << arg1 << " by " << arg2);
+    
+  private:
+                                    /**
+                                     * Prolongation matrix structures.
+                                     */
+    vector<SparseMatrixStruct > structures;
+                                    /**
+                                     * Prolongation matrices.
+                                     */
+    vector<SparseMatrix<float> > matrices;
+
+                                    /**
+                                     * Matrix generator.
+                                     * The arguments #nx# and #ny#
+                                     * are the numbers on the COARSE level.
+                                     * Builds a transfer matrix from
+                                     * fine to coarse (#vmult#).
+                                     */
+    void build_matrix(unsigned int nx, unsigned int ny,
+                     SparseMatrixStruct&, SparseMatrix<float>&);
+};

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.