]> https://gitweb.dealii.org/ - dealii.git/commitdiff
working on inconsistency
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 11 Nov 2009 13:23:09 +0000 (13:23 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 11 Nov 2009 13:23:09 +0000 (13:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@20089 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-39/step-39.cc

index 638037436f4b5cc1ae8284b57478fa5db8e01d0b..1ab4b2f563fe78bd63a93feec886ff2346e00597 100644 (file)
 #include <lac/compressed_sparsity_pattern.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
+#include <lac/precondition_block.h>
 
                                 // Include files for setting up the
                                 // mesh
 #include <grid/grid_generator.h>
 
-                                // Include files for FinitElement
+                                // Include files for FiniteElement
                                 // classes and DoFHandler.
+#include <fe/fe_q.h>
 #include <fe/fe_dgq.h>
 #include <dofs/dof_tools.h>
 #include <multigrid/mg_dof_handler.h>
 #include <numerics/mesh_worker.h>
 #include <numerics/mesh_worker_loop.h>
 
+                                // Support for multigrid methods
+#include <multigrid/mg_tools.h>
+#include <multigrid/multigrid.h>
+#include <multigrid/mg_matrix.h>
+#include <multigrid/mg_transfer.h>
+#include <multigrid/mg_coarse.h>
+#include <multigrid/mg_smoother.h>
+
                                 // Finally, we take our exact
-                                // solution from the library.
+                                // solution from the library as well
+                                // as quadrature and additional tools.
 #include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <numerics/vectors.h>
+#include <numerics/data_out.h>
 
 #include <iostream>
 #include <fstream>
                                 // as well.
 using namespace dealii;
 
+                                // This is the function we use to set
+                                // the boundary values and also the
+                                // exact solution we compare to.
+//Functions::LSingularityFunction exact_solution;
+Functions::Q1WedgeFunction<2> exact_solution;
 
                                 // @sect3{The local integrators}
 
@@ -177,16 +196,15 @@ void RHSIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceI
   const FEFaceValuesBase<dim>& fe = info.fe();
   Vector<double>& local_vector = info.R[0].block(0);
   
-  static Functions::LSingularityFunction lshaped;
   std::vector<double> boundary_values(fe.n_quadrature_points);
-  lshaped.value_list(fe.get_quadrature_points(), boundary_values);
+  exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
   
   const unsigned int deg = fe.get_fe().tensor_degree();
   const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
   
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
     for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-      local_vector(i) += (-fe.shape_value(i,k) * penalty * boundary_values[k]
+      local_vector(i) += (- fe.shape_value(i,k) * penalty * boundary_values[k]
                          + (fe.normal_vector(k) * fe.shape_grad(i,k)) * boundary_values[k])
                         * fe.JxW(k);
 }
@@ -210,8 +228,10 @@ class Step39
   private:
     void setup_system ();
     void assemble_matrix ();
+    void assemble_mg_matrix ();
     void assemble_right_hand_side ();
-    double estimate ();
+    void error ();
+    void estimate ();
     void solve ();
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
@@ -221,10 +241,16 @@ class Step39
     const FiniteElement<dim>& fe;
     MGDoFHandler<dim>      dof_handler;
 
-    SparsityPattern      sparsity_pattern;
+    SparsityPattern      sparsity;
     SparseMatrix<double> matrix;
     Vector<double>       solution;
     Vector<double>       right_hand_side;
+
+    MGLevelObject<SparsityPattern> mg_sparsity;
+    MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
+    MGLevelObject<SparseMatrix<double> > mg_matrix;
+    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_up;
+    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_down;
 };
 
 
@@ -234,7 +260,7 @@ Step39<dim>::Step39(const FiniteElement<dim>& fe)
                fe(fe),
                dof_handler(triangulation)
 {
-  GridGenerator::hyper_L(triangulation);
+  GridGenerator::hyper_cube(triangulation, -1, 1);
 }
 
 
@@ -248,12 +274,46 @@ Step39<dim>::setup_system()
   CompressedSparsityPattern c_sparsity(n_dofs);
   const DoFHandler<dim>& dof = dof_handler;
   DoFTools::make_flux_sparsity_pattern(dof, c_sparsity);
-  sparsity_pattern.copy_from(c_sparsity);
+  sparsity.copy_from(c_sparsity);
   
-  matrix.reinit(sparsity_pattern);
+  matrix.reinit(sparsity);
   
   solution.reinit(n_dofs);
   right_hand_side.reinit(n_dofs);
+
+  const unsigned int n_levels = triangulation.n_levels();
+  
+  mg_matrix.resize(0, n_levels-1);
+  mg_matrix.clear();
+  mg_matrix_dg_up.resize(0, n_levels-1);
+  mg_matrix_dg_up.clear();
+  mg_matrix_dg_down.resize(0, n_levels-1);
+  mg_matrix_dg_down.clear();
+
+  mg_sparsity.resize(0, n_levels-1);
+  mg_sparsity_dg_interface.resize(0, n_levels-1);
+  
+  for (unsigned int level=mg_sparsity.get_minlevel();
+       level<=mg_sparsity.get_maxlevel();++level)
+    {
+      CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level));
+      CompressedSparsityPattern ci_sparsity;
+      if (level>0)
+       ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level));
+      
+      MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level);
+      if (level>0)
+       MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level);
+      
+      mg_sparsity[level].copy_from(c_sparsity);
+      mg_matrix[level].reinit(mg_sparsity[level]);
+      if (level>0)
+       {
+         mg_sparsity_dg_interface[level].copy_from(ci_sparsity);
+         mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
+         mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
+       }
+    }
 }
 
 
@@ -276,6 +336,25 @@ Step39<dim>::assemble_matrix()
 }
 
 
+template <int dim>
+void
+Step39<dim>::assemble_mg_matrix()
+{
+  const MatrixIntegrator<dim> local;
+  MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> >, MatrixIntegrator<dim> >
+    integrator(local);
+  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+  integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+  UpdateFlags update_flags = update_values | update_gradients;
+  integrator.add_update_flags(update_flags, true, true, true, true);
+
+  integrator.initialize(mg_matrix);
+  MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
+  info_box.initialize(integrator, fe, mapping);
+  MeshWorker::integration_loop(dof_handler.begin(), dof_handler.end(), info_box, integrator);
+}
+
+
 template <int dim>
 void
 Step39<dim>::assemble_right_hand_side()
@@ -284,7 +363,7 @@ Step39<dim>::assemble_right_hand_side()
   MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::ResidualSimple<Vector<double> >, RHSIntegrator<dim> >
   integrator(local);
   const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
-  integrator.initialize_gauss_quadrature(1, n_gauss_points, n_gauss_points);
+  integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
   UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
   integrator.add_update_flags(update_flags, true, true, true, true);
 
@@ -295,6 +374,8 @@ Step39<dim>::assemble_right_hand_side()
   MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
   info_box.initialize(integrator, fe, mapping);
   MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator);
+  
+  right_hand_side *= -1.;
 }
 
 
@@ -304,7 +385,94 @@ Step39<dim>::solve()
 {
   SolverControl control(1000, 1.e-12);
   SolverCG<Vector<double> > cg(control);
-  cg.solve(matrix, solution, right_hand_side, PreconditionIdentity());
+
+  GrowingVectorMemory<Vector<double> > mem;
+  MGTransferPrebuilt<Vector<double> > mg_transfer;
+  mg_transfer.build_matrices(dof_handler);
+  FullMatrix<double> coarse_matrix;
+  coarse_matrix.copy_from (mg_matrix[0]);
+  MGCoarseGridHouseholder<double, Vector<double> > mg_coarse;
+  mg_coarse.initialize(coarse_matrix);
+  typedef PreconditionSOR<SparseMatrix<double> > RELAXATION;
+  MGSmootherRelaxation<SparseMatrix<double>, RELAXATION, Vector<double> >
+    mg_smoother(mem);
+  RELAXATION::AdditionalData smoother_data(1.);
+//  RELAXATION::AdditionalData smoother_data(fe.dofs_per_cell, 1.);
+  mg_smoother.initialize(mg_matrix, smoother_data);
+  
+                                  // Do two smoothing steps per level
+  mg_smoother.set_steps(2);
+                                  // Since the SOR method is not
+                                  // symmetric, but we use conjugate
+                                  // gradient iteration below, here
+                                  // is a trick to make the
+                                  // multilevel preconditioner a
+                                  // symmetric operator even for
+                                  // nonsymmetric smoothers.
+  mg_smoother.set_symmetric(true);
+  mg_smoother.set_variable(false);
+
+                                  // We must wrap our matrices in an
+                                  // object having the required
+                                  // multiplication functions.
+  MGMatrix<SparseMatrix<double>, Vector<double> > mgmatrix(&mg_matrix);
+  MGMatrix<SparseMatrix<double>, Vector<double> > mgdown(&mg_matrix_dg_down);
+  MGMatrix<SparseMatrix<double>, Vector<double> > mgup(&mg_matrix_dg_up);
+
+  
+                                  // Now, we are ready to set up the
+                                  // V-cycle operator and the
+                                  // multilevel preconditioner.
+  Multigrid<Vector<double> > mg(dof_handler, mgmatrix,
+                               mg_coarse, mg_transfer,
+                               mg_smoother, mg_smoother);
+  mg.set_edge_matrices(mgdown, mgup);
+  mg.set_debug(0);
+  mg_smoother.set_debug(0);
+  
+  PreconditionMG<dim, Vector<double>,
+    MGTransferPrebuilt<Vector<double> > >
+    preconditioner(dof_handler, mg, mg_transfer);
+  
+  cg.solve(matrix, solution, right_hand_side, preconditioner);
+}
+
+
+template <int dim>
+void
+Step39<dim>::error()
+{
+  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+2;
+  Vector<double> cell_errors(triangulation.n_active_cells());
+  
+  QGauss<dim> quadrature(n_gauss_points);
+  VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
+                                   cell_errors, quadrature, VectorTools::L2_norm);
+  deallog << "Error " << cell_errors.l2_norm() << std::endl;
+}
+
+
+template <int dim>
+void Step39<dim>::output_results (const unsigned int cycle) const
+{
+                                  // Output of the solution in
+                                  // gnuplot format.
+  std::string filename = "sol-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".gnuplot";
+  std::cout << "Writing solution to <" << filename << ">..."
+           << std::endl << std::endl;
+  std::ofstream gnuplot_output (filename.c_str());
+  
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "u");
+
+  data_out.build_patches ();
+  
+  data_out.write_gnuplot(gnuplot_output);
 }
 
 
@@ -329,15 +497,19 @@ Step39<dim>::run(unsigned int n_steps)
       deallog << std::endl;
 
       assemble_matrix();
+      assemble_mg_matrix();
       assemble_right_hand_side();
+      
       solve();
+      error();
+      output_results(s);
     }
 }
 
 
 int main()
 {
-  FE_DGQ<2> dgq1(1);
+  FE_Q<2> dgq1(1);
   Step39<2> test1(dgq1);
   test1.run(7);
 }

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.