]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Multigrid on meshes without hanging nodes
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 10 May 1999 12:57:34 +0000 (12:57 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 10 May 1999 12:57:34 +0000 (12:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@1307 0785d39b-7218-0410-832d-ea1e28bc413d

21 files changed:
deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/include/numerics/multigrid.h
deal.II/deal.II/include/numerics/multigrid.templates.h
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/deal.II/source/multigrid/mg_smoother.cc
deal.II/deal.II/source/multigrid/multigrid.cc
deal.II/deal.II/source/numerics/mg_smoother.cc
deal.II/deal.II/source/numerics/multigrid.cc
deal.II/lac/include/lac/mgbase.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/source/mgbase.cc
tests/deal.II/gradients.cc
tests/deal.II/helmholtz.h [new file with mode: 0644]
tests/deal.II/helmholtz1.th [new file with mode: 0644]
tests/deal.II/helmholtz1mg.th [new file with mode: 0644]
tests/deal.II/mg.cc [new file with mode: 0644]
tests/deal.II/mg2.cc [deleted file]
tests/deal.II/second_derivatives.cc
tests/lac/mg.cc

index 5f16cbf0e9d82346c0057f58213c478d8d6814ed..41c746197e3ba4e7faccf628f9b8b366333776a5 100644 (file)
@@ -318,6 +318,12 @@ class MGBase
                                      * Prolongation and restriction object.
                                      */
     SmartPointer<const MGTransferBase> transfer;
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
     
   protected:
 
@@ -352,7 +358,7 @@ class MGBase
                                      * Constructor, subject to change.
                                      */
     MGBase(const MGTransferBase& transfer,
-                 unsigned int maxlevel, unsigned int minlevel);
+                 unsigned int minlevel, unsigned int maxlevel);
                                     /**
                                      * Virtual destructor.
                                      */
index 916e2e7c6bc099123a7f99428c6b3d06143c1736..d9587dadfee635d179fb3f52c54777bf3503b440 100644 (file)
@@ -28,7 +28,9 @@ class MG
   public MGBase
 {
   public:
-    MG(const MGDoFHandler<dim>&, const MGTransferBase& transfer);
+    MG(const MGDoFHandler<dim>&,
+       const MGMatrix<SparseMatrix<float> >&,
+       const MGTransferBase& transfer);
     
                                     /** Transfer from dVector to
                                      * MGVector.
@@ -40,8 +42,7 @@ class MG
                                      * MGVector.
                                      */
     template<typename number>
-    void copy_to_mg(vector<Vector<float> >& dst,
-                   const Vector<number>& src) const;
+    void copy_to_mg(const Vector<number>& src);
 
                                     /**
                                      * Transfer from multi-grid vector to
@@ -54,23 +55,15 @@ class MG
                                      * #src# are zero.
                                      */
     template<typename number>
-    void copy_from_mg(Vector<number>& dst,
-                     const vector<Vector<float> >& src) const;
+    void copy_from_mg(Vector<number>& dst) const;
 
                                     /**
-                                     * Access to matrix structure.
+                                     * Residual on a level.
                                      */
-    SparseMatrixStruct& get_sparsity_pattern(unsigned int level);
-
-                                    /**
-                                     * Read-only access to matrix structure.
-                                     */
-    const SparseMatrixStruct& get_sparsity_pattern(unsigned int level) const;
-
-                                    /**
-                                     * Access to level matrices.
-                                     */
-    SparseMatrix<float>& get_matrix(unsigned int level);
+    virtual void level_residual(unsigned int,
+                               Vector<float> &,
+                               const Vector<float> &,
+                               const Vector<float> &);
 
                                     /**
                                      * Read-only access to level matrices.
@@ -83,21 +76,13 @@ class MG
                                      */
     SmartPointer<const MGDoFHandler<dim> > dofs;
 
-                                    /**
-                                     * Matrix structures for each level.
-                                     * These are generated by the constructor
-                                     * of #MG# and can be accessed
-                                     * later.
-                                     */
-    vector<SparseMatrixStruct> matrix_structures;
-
                                     /**
                                      * Matrices for each level.
                                      * The matrices are prepared by
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    vector<SparseMatrix<float> > matrices;
+    SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
 };
 
 
@@ -175,16 +160,6 @@ class MGTransferPrebuilt : public MGTransferBase
 };
 
 
-template<int dim>
-inline
-SparseMatrix<float>&
-MG<dim>::get_matrix(unsigned int level)
-{
-  Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
-  
-  return matrices[level];
-}
-
 template<int dim>
 inline
 const SparseMatrix<float>&
@@ -192,7 +167,7 @@ MG<dim>::get_matrix(unsigned int level) const
 {
   Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
   
-  return matrices[level];
+  return (*matrices)[level];
 }
 
 
index d290cd72238bbeac47163282ebc8d0cede109ac0..9e170e3d5b2553474f24f346d41d49daca682d4d 100644 (file)
@@ -2,42 +2,72 @@
 
 #include <numerics/multigrid.h>
 
-/*
 template <int dim>
-MultiGrid::
-*/
-
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs, const MGTransferBase& transfer)
+MG<dim>::MG(const MGDoFHandler<dim>& dofs,
+           const MGMatrix<SparseMatrix<float> >& matrices,
+           const MGTransferBase& transfer)
                :
                MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
                dofs(&dofs),
-               matrix_structures(maxlevel-minlevel),
-               matrices(maxlevel-minlevel)
+               matrices(&matrices)
 {
-  for(unsigned l = minlevel; l < maxlevel; ++l)
+  for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
-      dofs.make_sparsity_pattern(l, matrix_structures[l-minlevel]);
-      matrices[l-minlevel].reinit(matrix_structures[l-minlevel]);
+      s[level].reinit(matrices[level].m());
+      d[level].reinit(matrices[level].m());
     }
 }
 
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
 void
-MG<dim>::copy_to_mg(vector<Vector<float> >& dst,
-                          const Vector<number>& src) const
+MG<dim>::copy_to_mg(const Vector<number>& src)
 {
+  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  
+  vector<int> index(fe_dofs);
+  vector<int> mgindex(fe_dofs);
 
+  DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
+                                    ; c != dofs->end() ; ++c, ++dc)
+    {
+      dc->get_dof_indices(index);
+      c->get_mg_dof_indices(mgindex);
+      for (unsigned int i=0;i<fe_dofs;++i)
+       d[maxlevel](mgindex[i]) = src(index[i]);
+    } 
 }
 
 template <int dim> template <typename number>
 void
-MG<dim>::copy_from_mg(Vector<number>& dst,
-                            const vector<Vector<float> >& src) const
+MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
-// active iterator
+  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+  vector<int> index(fe_dofs);
+  vector<int> mgindex(fe_dofs);
+
+  DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
+                                    ; c != dofs->end() ; ++c, ++dc)
+    {
+      dc->get_dof_indices(index);
+      c->get_mg_dof_indices(mgindex);
+      for (unsigned int i=0;i<fe_dofs;++i)
+       dst(index[i]) = s[maxlevel](mgindex[i]);
+    }
 }
 
+template <int dim>
+void
+MG<dim>::level_residual(unsigned int l,
+                       Vector<float> &r,
+                       const Vector<float> &u,
+                       const Vector<float> &b)
+{
+  (*matrices)[l].residual(r,u,b);
+}
 
 /*
 template <int dim>
@@ -46,7 +76,5 @@ template <int dim>
 MG::
 template <int dim>
 MG::
-template <int dim>
-MG::
 
 */
index 916e2e7c6bc099123a7f99428c6b3d06143c1736..d9587dadfee635d179fb3f52c54777bf3503b440 100644 (file)
@@ -28,7 +28,9 @@ class MG
   public MGBase
 {
   public:
-    MG(const MGDoFHandler<dim>&, const MGTransferBase& transfer);
+    MG(const MGDoFHandler<dim>&,
+       const MGMatrix<SparseMatrix<float> >&,
+       const MGTransferBase& transfer);
     
                                     /** Transfer from dVector to
                                      * MGVector.
@@ -40,8 +42,7 @@ class MG
                                      * MGVector.
                                      */
     template<typename number>
-    void copy_to_mg(vector<Vector<float> >& dst,
-                   const Vector<number>& src) const;
+    void copy_to_mg(const Vector<number>& src);
 
                                     /**
                                      * Transfer from multi-grid vector to
@@ -54,23 +55,15 @@ class MG
                                      * #src# are zero.
                                      */
     template<typename number>
-    void copy_from_mg(Vector<number>& dst,
-                     const vector<Vector<float> >& src) const;
+    void copy_from_mg(Vector<number>& dst) const;
 
                                     /**
-                                     * Access to matrix structure.
+                                     * Residual on a level.
                                      */
-    SparseMatrixStruct& get_sparsity_pattern(unsigned int level);
-
-                                    /**
-                                     * Read-only access to matrix structure.
-                                     */
-    const SparseMatrixStruct& get_sparsity_pattern(unsigned int level) const;
-
-                                    /**
-                                     * Access to level matrices.
-                                     */
-    SparseMatrix<float>& get_matrix(unsigned int level);
+    virtual void level_residual(unsigned int,
+                               Vector<float> &,
+                               const Vector<float> &,
+                               const Vector<float> &);
 
                                     /**
                                      * Read-only access to level matrices.
@@ -83,21 +76,13 @@ class MG
                                      */
     SmartPointer<const MGDoFHandler<dim> > dofs;
 
-                                    /**
-                                     * Matrix structures for each level.
-                                     * These are generated by the constructor
-                                     * of #MG# and can be accessed
-                                     * later.
-                                     */
-    vector<SparseMatrixStruct> matrix_structures;
-
                                     /**
                                      * Matrices for each level.
                                      * The matrices are prepared by
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    vector<SparseMatrix<float> > matrices;
+    SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
 };
 
 
@@ -175,16 +160,6 @@ class MGTransferPrebuilt : public MGTransferBase
 };
 
 
-template<int dim>
-inline
-SparseMatrix<float>&
-MG<dim>::get_matrix(unsigned int level)
-{
-  Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
-  
-  return matrices[level];
-}
-
 template<int dim>
 inline
 const SparseMatrix<float>&
@@ -192,7 +167,7 @@ MG<dim>::get_matrix(unsigned int level) const
 {
   Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
   
-  return matrices[level];
+  return (*matrices)[level];
 }
 
 
index d290cd72238bbeac47163282ebc8d0cede109ac0..9e170e3d5b2553474f24f346d41d49daca682d4d 100644 (file)
@@ -2,42 +2,72 @@
 
 #include <numerics/multigrid.h>
 
-/*
 template <int dim>
-MultiGrid::
-*/
-
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs, const MGTransferBase& transfer)
+MG<dim>::MG(const MGDoFHandler<dim>& dofs,
+           const MGMatrix<SparseMatrix<float> >& matrices,
+           const MGTransferBase& transfer)
                :
                MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
                dofs(&dofs),
-               matrix_structures(maxlevel-minlevel),
-               matrices(maxlevel-minlevel)
+               matrices(&matrices)
 {
-  for(unsigned l = minlevel; l < maxlevel; ++l)
+  for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
-      dofs.make_sparsity_pattern(l, matrix_structures[l-minlevel]);
-      matrices[l-minlevel].reinit(matrix_structures[l-minlevel]);
+      s[level].reinit(matrices[level].m());
+      d[level].reinit(matrices[level].m());
     }
 }
 
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
 void
-MG<dim>::copy_to_mg(vector<Vector<float> >& dst,
-                          const Vector<number>& src) const
+MG<dim>::copy_to_mg(const Vector<number>& src)
 {
+  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  
+  vector<int> index(fe_dofs);
+  vector<int> mgindex(fe_dofs);
 
+  DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
+                                    ; c != dofs->end() ; ++c, ++dc)
+    {
+      dc->get_dof_indices(index);
+      c->get_mg_dof_indices(mgindex);
+      for (unsigned int i=0;i<fe_dofs;++i)
+       d[maxlevel](mgindex[i]) = src(index[i]);
+    } 
 }
 
 template <int dim> template <typename number>
 void
-MG<dim>::copy_from_mg(Vector<number>& dst,
-                            const vector<Vector<float> >& src) const
+MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
-// active iterator
+  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+  vector<int> index(fe_dofs);
+  vector<int> mgindex(fe_dofs);
+
+  DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
+                                    ; c != dofs->end() ; ++c, ++dc)
+    {
+      dc->get_dof_indices(index);
+      c->get_mg_dof_indices(mgindex);
+      for (unsigned int i=0;i<fe_dofs;++i)
+       dst(index[i]) = s[maxlevel](mgindex[i]);
+    }
 }
 
+template <int dim>
+void
+MG<dim>::level_residual(unsigned int l,
+                       Vector<float> &r,
+                       const Vector<float> &u,
+                       const Vector<float> &b)
+{
+  (*matrices)[l].residual(r,u,b);
+}
 
 /*
 template <int dim>
@@ -46,7 +76,5 @@ template <int dim>
 MG::
 template <int dim>
 MG::
-template <int dim>
-MG::
 
 */
index b08359efb40eba1252c8022d3f9dab6e10d1fdea..e6340762e8d44d4eb285bd9bfb3952abbd780e4f 100644 (file)
@@ -26,13 +26,14 @@ MGBase::~MGBase ()
 
 
 MGBase::MGBase(const MGTransferBase& transfer,
-                            unsigned maxlevel, unsigned minlevel)
+                            unsigned minlevel, unsigned maxlevel)
                :
                maxlevel(maxlevel), minlevel(minlevel),
                d(minlevel,maxlevel),
                s(minlevel,maxlevel),
                transfer(&transfer)
 {
+  Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
 }
 
 void
index 71b80cb4b3c4af2b7454baf00a08a7b011c6eea2..23bea7acf92eb1e14eb0ed1197fe88f1e87a8427 100644 (file)
@@ -90,11 +90,6 @@ MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
 
 //////////////////////////////////////////////////////////////////////
 
-MGSmootherBase::~MGSmootherBase () 
-{};
-
-//////////////////////////////////////////////////////////////////////
-
 void
 MGSmoother::set_zero_interior_boundary (const unsigned int  level,
                                        Vector<float>      &u) const
index b5db616ddeaf3f837de3f5ba4cf0e21647cde99f..429f929615908d12effdc73a45f973ac2543c2ce 100644 (file)
@@ -137,8 +137,8 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
                                     Vector<float>       &dst,
                                     const Vector<float> &src) const 
 {
-  Assert ((to_level >= 1) && (to_level<prolongation_matrices.size()-1),
-         ExcIndexRange (to_level, 1, prolongation_matrices.size()-1));
+  Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+         ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
 
   prolongation_matrices[to_level-1].vmult (dst, src);
 };
@@ -149,8 +149,8 @@ void MGTransferPrebuilt::restrict (const unsigned int   from_level,
                                   Vector<float>       &dst,
                                   const Vector<float> &src) const 
 {
-  Assert ((from_level >= 1) && (from_level<prolongation_matrices.size()-1),
-         ExcIndexRange (from_level, 1, prolongation_matrices.size()-1));
+  Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+         ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
 
   prolongation_matrices[from_level-1].Tvmult (dst, src);
 };
index 71b80cb4b3c4af2b7454baf00a08a7b011c6eea2..23bea7acf92eb1e14eb0ed1197fe88f1e87a8427 100644 (file)
@@ -90,11 +90,6 @@ MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
 
 //////////////////////////////////////////////////////////////////////
 
-MGSmootherBase::~MGSmootherBase () 
-{};
-
-//////////////////////////////////////////////////////////////////////
-
 void
 MGSmoother::set_zero_interior_boundary (const unsigned int  level,
                                        Vector<float>      &u) const
index b5db616ddeaf3f837de3f5ba4cf0e21647cde99f..429f929615908d12effdc73a45f973ac2543c2ce 100644 (file)
@@ -137,8 +137,8 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
                                     Vector<float>       &dst,
                                     const Vector<float> &src) const 
 {
-  Assert ((to_level >= 1) && (to_level<prolongation_matrices.size()-1),
-         ExcIndexRange (to_level, 1, prolongation_matrices.size()-1));
+  Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+         ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
 
   prolongation_matrices[to_level-1].vmult (dst, src);
 };
@@ -149,8 +149,8 @@ void MGTransferPrebuilt::restrict (const unsigned int   from_level,
                                   Vector<float>       &dst,
                                   const Vector<float> &src) const 
 {
-  Assert ((from_level >= 1) && (from_level<prolongation_matrices.size()-1),
-         ExcIndexRange (from_level, 1, prolongation_matrices.size()-1));
+  Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+         ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
 
   prolongation_matrices[from_level-1].Tvmult (dst, src);
 };
index 3e9396f8f0ee5e467308e41d379dbf032d293346..cc71bc98f4dcadfb6ea3394de9124fef0720f140 100644 (file)
@@ -318,6 +318,12 @@ class MGBase
                                      * Prolongation and restriction object.
                                      */
     SmartPointer<const MGTransferBase> transfer;
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
     
   protected:
 
@@ -352,7 +358,7 @@ class MGBase
                                      * Constructor, subject to change.
                                      */
     MGBase(const MGTransferBase& transfer,
-                 unsigned int maxlevel, unsigned int minlevel);
+                 unsigned int minlevel, unsigned int maxlevel);
                                     /**
                                      * Virtual destructor.
                                      */
index a2342846fdf96a88b663fc4cd548c0c8767e5547..5218fd84b0cd6faa8b004b994bb1733165a6500a 100644 (file)
@@ -489,7 +489,9 @@ SparseMatrix<number>::SSOR (Vector<somenumber>& dst, const number om) const
 
 
 template <typename number>
-const SparseMatrixStruct & SparseMatrix<number>::get_sparsity_pattern () const {
+const SparseMatrixStruct &
+SparseMatrix<number>::get_sparsity_pattern () const
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   return *cols;
 };
index b08359efb40eba1252c8022d3f9dab6e10d1fdea..e6340762e8d44d4eb285bd9bfb3952abbd780e4f 100644 (file)
@@ -26,13 +26,14 @@ MGBase::~MGBase ()
 
 
 MGBase::MGBase(const MGTransferBase& transfer,
-                            unsigned maxlevel, unsigned minlevel)
+                            unsigned minlevel, unsigned maxlevel)
                :
                maxlevel(maxlevel), minlevel(minlevel),
                d(minlevel,maxlevel),
                s(minlevel,maxlevel),
                transfer(&transfer)
 {
+  Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
 }
 
 void
index 3963071d7387fad26333dd7dfb2419beb72f8cec..1b000fc53e40e8bdedaf265be0598c780d3d1fe6 100644 (file)
@@ -36,7 +36,7 @@ int main () {
   dof.distribute_dofs(fe);
 
   QTrapez<2> q;
-  FEValues<2> fevalues(fe,q,update_gradients,tria.get_boundary());
+  FEValues<2> fevalues(fe,q,update_gradients);
   fevalues.reinit (dof.begin_active());
   
   
diff --git a/tests/deal.II/helmholtz.h b/tests/deal.II/helmholtz.h
new file mode 100644 (file)
index 0000000..f511560
--- /dev/null
@@ -0,0 +1,26 @@
+// $Id$
+
+#include <lac/mgbase.h>
+
+/**
+ * Generator for system matrix and right hand side.
+ */
+class Helmholtz
+{
+  public:
+    template<int dim, class MATRIX, class VECTOR>
+    void build_all(MATRIX& A,
+                  VECTOR& f,
+                  const DoFHandler<dim>& dof,
+                  const Quadrature<dim>& quadrature,
+                  const Function<dim>& rhs);
+
+    template<int dim, class MATRIX>
+    void build_mgmatrix(MGMatrix<MATRIX>& A,
+                  const MGDoFHandler<dim>& dof,
+                  const Quadrature<dim>& quadrature);
+};
+
+#include "helmholtz1.th"
+#include "helmholtz1mg.th"
+
diff --git a/tests/deal.II/helmholtz1.th b/tests/deal.II/helmholtz1.th
new file mode 100644 (file)
index 0000000..5abed88
--- /dev/null
@@ -0,0 +1,70 @@
+// $Id$
+
+template<int dim, class MATRIX, class VECTOR>
+void
+Helmholtz::build_all(MATRIX& A,
+                    VECTOR& f,
+                    const DoFHandler<dim>& dof,
+                    const Quadrature<dim>& quadrature,
+                    const Function<dim>& rhs)
+{
+  const unsigned int fe_dofs = dof.get_fe().total_dofs;
+  FEValues<dim> fe(dof.get_fe(), quadrature,
+                  UpdateFlags(update_gradients | update_JxW_values |
+                              update_q_points));
+  
+  Vector<double> elvec(fe_dofs);
+  FullMatrix<double> elmat(fe_dofs);
+  vector<int> indices(fe_dofs);
+
+////////////////////////////////////////////////////
+// Loop for computing the element matrices
+////////////////////////////////////////////////////
+    
+  for (DoFHandler<dim>::active_cell_iterator c = dof.begin_active()
+                                       ; c != dof.end() ; ++c)
+    {
+      elmat.clear();
+      elvec.clear();
+      fe.reinit(c);
+
+                                      // Quadrature loop      
+
+      for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
+       {
+         Point<dim> loc = fe.quadrature_point(k);
+
+                                          // Test function loop
+         for (unsigned i=0;i<fe_dofs;++i)
+           {
+             const Point<dim> dv = fe.shape_grad(i,k);
+             double v = fe.shape_value(i,k);
+             
+             elvec(i) += fe.JxW(k) *
+                       rhs(loc) * v;
+             
+                                              //Trial function loop
+             for (unsigned j=0;j<fe_dofs;++j)
+               {
+                 const Point<dim> du = fe.shape_grad(j,k);
+                 double u = fe.shape_value(j,k);
+                 
+                 elmat(i,j) += fe.JxW(k) *
+                               (1.*u*v + du*dv)
+                               ;
+               }
+           }
+       }
+
+      c->get_dof_indices(indices);
+      for (unsigned i=0;i<fe_dofs;++i)
+       {
+         f(indices[i]) += elvec(i);
+         
+         for (unsigned j=0;j<fe_dofs;++j)
+           {
+             A.add(indices[i], indices[j], elmat(i,j));
+           }
+       }
+    }
+}
diff --git a/tests/deal.II/helmholtz1mg.th b/tests/deal.II/helmholtz1mg.th
new file mode 100644 (file)
index 0000000..450e0a4
--- /dev/null
@@ -0,0 +1,66 @@
+// $Id$
+
+template<int dim, class MATRIX>
+void
+Helmholtz::build_mgmatrix(MGMatrix<MATRIX>& A,
+                         const MGDoFHandler<dim>& dof,
+                         const Quadrature<dim>& quadrature)
+{
+  const unsigned int fe_dofs = dof.get_fe().total_dofs;
+  FEValues<dim> fe(dof.get_fe(), quadrature,
+                  UpdateFlags(update_gradients | update_JxW_values |
+                              update_q_points));
+  
+  Vector<double> elvec(fe_dofs);
+  FullMatrix<double> elmat(fe_dofs);
+  vector<int> indices(fe_dofs);
+
+////////////////////////////////////////////////////
+// Loop for computing the element matrices
+////////////////////////////////////////////////////
+    
+    DoFHandler<dim>::cell_iterator dc = dof.DoFHandler<dim>::begin();
+    for (MGDoFHandler<dim>::cell_iterator c = dof.begin()
+                                    ; c != dof.end() ; ++c, ++dc)
+    {
+      elmat.clear();
+      elvec.clear();
+      fe.reinit(dc);
+
+                                      // Quadrature loop      
+
+      for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
+       {
+//       Point<dim> loc = fe.quadrature_point(k);
+
+                                          // Test function loop
+         for (unsigned i=0;i<fe_dofs;++i)
+           {
+             const Point<dim> dv = fe.shape_grad(i,k);
+             double v = fe.shape_value(i,k);
+             
+                                              //Trial function loop
+             for (unsigned j=0;j<fe_dofs;++j)
+               {
+                 const Point<dim> du = fe.shape_grad(j,k);
+                 double u = fe.shape_value(j,k);
+                 
+                 elmat(i,j) += fe.JxW(k) *
+                               (1.*u*v + du*dv)
+                               ;
+               }
+           }
+       }
+
+      const unsigned int level = c->level();
+      
+      c->get_mg_dof_indices(indices);
+      for (unsigned i=0;i<fe_dofs;++i)
+       {
+         for (unsigned j=0;j<fe_dofs;++j)
+           {
+             A[level].add(indices[i], indices[j], elmat(i,j));
+           }
+       }
+    }
+}
diff --git a/tests/deal.II/mg.cc b/tests/deal.II/mg.cc
new file mode 100644 (file)
index 0000000..29aa26d
--- /dev/null
@@ -0,0 +1,160 @@
+// $Id$
+
+// deal_II_libraries.g=-ldeal_II_2d.g
+// deal_II_libraries=-ldeal_II_2d
+
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <lac/sparsematrix.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_richardson.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_accessor.h>
+#include <grid/mg_dof_accessor.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_values.h>
+#include <numerics/multigrid.h>
+#include <numerics/multigrid.templates.h>
+#include <numerics/mg_smoother.h>
+
+#include "helmholtz.h"
+
+template<int dim>
+class RHSFunction
+  :
+  public Function<dim>
+{
+  public:
+    virtual double operator() (const Point<dim>&) const;
+//    virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+};
+
+class MGSmootherLAC
+  :
+  public 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;
+    
+};
+
+main(int argc, const char** argv)
+{
+  unsigned int refine;
+  if (argc<=1)
+    {
+      cerr << "Usage:" << argv[0] << " <refinement>" << endl;
+      refine = 0;
+    }
+  else
+    refine = atoi(argv[1]);
+  
+
+  Helmholtz equation;
+  RHSFunction<2> rhs;
+  QGauss3<2> quadrature;
+  
+  FELinear<2> fe;
+  Triangulation<2> tr;
+  MGDoFHandler<2> mgdof(&tr);
+  DoFHandler<2>& dof(mgdof);
+  
+  GridGenerator::hyper_cube(tr,-1.,1.);
+  tr.refine_global(refine);
+  dof.distribute_dofs(fe);
+
+  const unsigned int size = dof.n_dofs();
+  
+  SparseMatrixStruct structure(size, dof.max_couplings_between_dofs());
+  dof.make_sparsity_pattern(structure);
+  structure.compress();
+  
+  SparseMatrix<double> A(structure);
+  Vector<double> f(size);
+
+  equation.build_all(A,f,dof, quadrature, rhs);
+
+  Vector<double> u;
+  u.reinit(f);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverControl control(100, 1.e-15);
+  SolverCG<SparseMatrix<double>, Vector<double> >solver(control, mem);
+  
+  PreconditionIdentity<Vector<double> > precondition;
+  solver.solve(A,u,f,precondition);
+  
+  u = 0.;
+  vector<SparseMatrixStruct> mgstruct(tr.n_levels());
+  MGMatrix<SparseMatrix<float> > mgA(0,tr.n_levels()-1);
+  for (unsigned int i=0;i<tr.n_levels();++i)
+    {
+      mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
+                        mgdof.max_couplings_between_dofs());
+      mgdof.make_sparsity_pattern(i, mgstruct[i]);
+      mgstruct[i].compress();
+
+      mgA[i].reinit(mgstruct[i]);
+    }
+  equation.build_mgmatrix(mgA, mgdof, quadrature);
+  
+  SolverControl cgcontrol(10,0., false, false);
+  
+  PrimitiveVectorMemory<Vector<float> > cgmem;
+  SolverCG<SparseMatrix<float>, Vector<float> > cgsolver(cgcontrol, cgmem);
+  PreconditionIdentity<Vector<float> > cgprec;
+  MGCoarseGridLACIteration<SolverCG<SparseMatrix<float>, Vector<float> >,
+    SparseMatrix<float>, PreconditionIdentity<Vector<float> > >
+    coarse(cgsolver, mgA[0], cgprec);
+
+  MGSmootherLAC smoother(mgA);
+  MGTransferPrebuilt transfer;
+  transfer.build_matrices(mgdof);
+  
+  deallog << "Levels: " << tr.n_levels() << endl;
+  
+  MG<2> multigrid(mgdof, mgA, transfer);
+  PreconditionMG<MG<2>, Vector<double> >
+    mgprecondition(multigrid, smoother, smoother, coarse);
+  
+  solver.solve(A, u, f, mgprecondition);
+  
+  cerr << "OK" << endl;
+}
+
+template<int dim>
+double
+RHSFunction<dim>::operator() (const Point<dim>&) const
+{
+  return 1.;
+}
+
+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(1,1.e-300,false,false);
+  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.);
+
+  rich.solve((*matrices)[level], u, rhs, prec);
+}
diff --git a/tests/deal.II/mg2.cc b/tests/deal.II/mg2.cc
deleted file mode 100644 (file)
index 44dd281..0000000
+++ /dev/null
@@ -1,23 +0,0 @@
-// $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 <grid/grid_generator.h>
-#include <fe/fe_lib.lagrange.h>
-
-main()
-{
-  Triangulation<2> tr;
-  MGDoFHandler<2> dof(&tr);
-  FELinear<2> fe;
-  
-  GridGenerator::hyper_cube(tr,-1.,1.);
-  tr.refine_global(3);
-  dof.distribute_dofs(fe);
-}
-
index 2ad62be46c59fdf8cd29874bdfc63d360ed3e542..0c4d5064cb95e9fafd9082b0638426b9f02913dd 100644 (file)
@@ -28,7 +28,7 @@ int main () {
 
   StraightBoundary<2> b;
   QTrapez<2> q;
-  FEValues<2> fevalues(fe,q,update_second_derivatives,b);
+  FEValues<2> fevalues(fe,q,update_second_derivatives);
   
   
   Vector<double> val(4);
index a8df6709d660cce5e03d81f70658bba4a98f3315..c667839d0b4aa1865a64debb7bae61f6eb4f4f07 100644 (file)
@@ -160,7 +160,7 @@ main()
 FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
           FDMGTransfer& transfer)
                :
-               MGBase(transfer, maxlevel, 0),
+               MGBase(transfer, 0, maxlevel),
                matrices(&matrices)
 {
   for (unsigned int level = minlevel; level<=maxlevel ; ++level)

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.