]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New test, which L_2 projects a polynomial function
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 29 Jul 2006 23:43:28 +0000 (23:43 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 29 Jul 2006 23:43:28 +0000 (23:43 +0000)
to a hp FE discretisation and afterwards evaluates
the L_2 norm of the error. The L_2 norm of the error
should be in the range of the FP accuracy.

git-svn-id: https://svn.dealii.org/trunk@13513 0785d39b-7218-0410-832d-ea1e28bc413d

tests/hp/hp_hanging_nodes_02.cc [new file with mode: 0644]

diff --git a/tests/hp/hp_hanging_nodes_02.cc b/tests/hp/hp_hanging_nodes_02.cc
new file mode 100644 (file)
index 0000000..e41eff5
--- /dev/null
@@ -0,0 +1,417 @@
+//----------------------------  step-3.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 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.
+//
+//----------------------------  step-3.cc  ---------------------------
+
+
+// modified step-3 to be a L_2 projection. This function
+// L_2-projects a bilinear, biquadratic or bicubic function
+// onto a square shaped domain which is filled with
+// hp elements. Afterwards the L_2 error is computed, which
+// should be in the range of the FP approximation order.
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <grid/grid_generator.h>
+
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+
+#include <fe/fe_q.h>
+
+#include <dofs/dof_tools.h>
+
+#include <fe/hp_fe_values.h>
+#include <base/quadrature_lib.h>
+
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+
+std::ofstream logfile("hp_hanging_nodes_02/output");
+
+
+template <int dim>
+class Linear : public Function<dim> 
+{
+  public:
+    Linear () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;    
+};
+
+
+
+template <int dim>
+double Linear<dim>::value (const Point<dim> &p,
+                          const unsigned int) const 
+{
+  double res = p(0);
+  for (unsigned int d = 1; d < dim; ++d)
+    res *= p(d);
+  
+  return res;  
+}
+
+
+
+template <int dim>
+class Quadratic : public Function<dim> 
+{
+  public:
+    Quadratic () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;    
+};
+
+
+
+template <int dim>
+double Quadratic<dim>::value (const Point<dim> &p,
+                          const unsigned int) const 
+{
+  double res = p(0) * p(0);
+  for (unsigned int d = 1; d < dim; ++d)
+    res *= p(d) * p(d);
+  
+  return res;  
+}
+
+
+template <int dim>
+class Cubic : public Function<dim> 
+{
+  public:
+    Cubic () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;    
+};
+
+
+
+template <int dim>
+double Cubic<dim>::value (const Point<dim> &p,
+                          const unsigned int) const 
+{
+  double res = p(0) * p(0) * p(0);
+  for (unsigned int d = 1; d < dim; ++d)
+    res *= p(d) * p(d) * p(d);
+  
+  return res;  
+}
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+
+    void run (const Function<dim> &f_test,
+             bool random,
+             unsigned int *indx);
+    
+  private:
+    void make_grid_and_dofs (const bool random_p);
+    void assemble_system (const Function<dim> &f_test);
+    void solve ();
+    void eval_error (const Function<dim> &f_test);
+    void output_results () const;
+
+    Triangulation<dim>     triangulation;
+    hp::FECollection<dim>              fe;
+    hp::DoFHandler<dim>        dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+                                    // Although we do not have h-refinement,
+                                    // hanging nodes will inevitably appear
+                                    // due to different polynomial degrees.
+    ConstraintMatrix     hanging_node_constraints;
+    
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{}
+
+
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs (const bool random_p)
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (3);
+  deallog << "Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl;
+  deallog << "Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl;
+
+                                  // Now to the p-Method. Assign
+                                  // random active_fe_indices to the
+                                  // different cells.
+  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active (),
+                                                    endc = dof_handler.end ();
+  if (random_p)
+    {      
+      for (; cell != endc; ++cell)
+       {      
+         cell->set_active_fe_index ((int)(4.0 * (double) random () / (double) RAND_MAX));
+       }
+    }
+  else
+    {      
+      unsigned int cell_no = 0;  
+      for (; cell != endc; ++cell)
+       {
+         if (cell_no >= triangulation.n_active_cells () / 2)
+           cell->set_active_fe_index (1);
+         else
+           cell->set_active_fe_index (0);
+       }  
+    }
+
+  
+  dof_handler.distribute_dofs (fe);
+  deallog << "Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+                                  // Create sparsity pattern.  
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+                                  // Create constraints which stem from
+                                  // the different polynomial degrees on
+                                  // the different elements.
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          hanging_node_constraints);
+
+  hanging_node_constraints.print (deallog.get_file_stream ());
+
+  hanging_node_constraints.close ();
+  hanging_node_constraints.condense (sparsity_pattern);
+  
+  sparsity_pattern.compress();
+  system_matrix.reinit (sparsity_pattern);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system (const Function<dim> &f_test) 
+{
+  hp::QCollection<dim>  quadrature_formula(QGauss<dim>(6));
+  hp::FEValues<dim> x_fe_values (fe, quadrature_formula, 
+                        update_values | update_gradients |
+                              update_JxW_values | update_q_points);
+  
+  const unsigned int   max_dofs_per_cell = fe.max_dofs_per_cell ();
+  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (max_dofs_per_cell, max_dofs_per_cell);
+  Vector<double>       cell_rhs (max_dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (max_dofs_per_cell);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      x_fe_values.reinit (cell);
+
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+      
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      const unsigned int dofs_per_cell = cell->get_fe ().dofs_per_cell;
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+           cell_matrix(i,j) += (fe_values.shape_value (i, q_point) *
+                                fe_values.shape_value (j, q_point) *
+                                fe_values.JxW (q_point));
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         {         
+           Point<dim> p_q = fe_values.quadrature_point (q_point);
+           double f = f_test.value (p_q);          
+           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                           f *
+                           fe_values.JxW (q_point));
+         }      
+           
+      local_dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (local_dof_indices);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add (local_dof_indices[i],
+                            local_dof_indices[j],
+                            cell_matrix(i,j));
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       system_rhs(local_dof_indices[i]) += cell_rhs(i);
+    }
+
+                                  // Include hanging nodes.
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  cg.solve (system_matrix, solution, system_rhs,
+           PreconditionIdentity());
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::eval_error (const Function<dim> &f_test) 
+{
+  hp::QCollection<dim>  quadrature(QGauss<dim>(6));
+  Vector<double> cellwise_errors (triangulation.n_active_cells());
+  VectorTools::integrate_difference (dof_handler, solution, f_test,
+                                     cellwise_errors, quadrature,
+                                     VectorTools::L2_norm);
+  const double p_l2_error = cellwise_errors.l2_norm();
+
+  deallog << "L2_Error : " << p_l2_error << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results () const
+{
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+
+  data_out.write_gnuplot (deallog.get_file_stream());
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run (const Function<dim> &f_test,
+                              bool random,
+                              unsigned int *indx) 
+{
+  FE_Q<dim> fe_1 (indx[0]),
+    fe_2 (indx[1]),
+    fe_3 (indx[2]),
+    fe_4 (indx[3]);
+
+  fe.push_back (fe_1);
+  fe.push_back (fe_2);
+  fe.push_back (fe_3);
+  fe.push_back (fe_4);
+  
+  make_grid_and_dofs (random);
+  assemble_system (f_test);
+  solve ();
+  eval_error (f_test);
+  output_results ();
+}
+
+
+template <int dim>
+void run_test (const Function<dim> &f_test,
+              unsigned int *indx)
+{
+  LaplaceProblem<dim> laplace_problem_1;
+  laplace_problem_1.run (f_test, true, indx);
+
+  LaplaceProblem<dim> laplace_problem_2;
+  laplace_problem_2.run (f_test, false, indx);
+}
+
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  unsigned int index[] = 
+    {
+         1,2,3,4,5,6,7
+    };
+  
+  
+  Linear<2> test2d_1;
+  Quadratic<2> test2d_2;
+  Cubic<2> test2d_3;
+
+  Linear<3> test3d_1;
+  Quadratic<3> test3d_2;
+  Cubic<3> test3d_3;
+
+  deallog << "Testing Order 1" << std::endl;
+  run_test<2> (test2d_1, &(index[0]));
+  run_test<3> (test3d_1, &(index[0]));
+
+  deallog << "Testing Order 2" << std::endl;
+  run_test<2> (test2d_2, &(index[1]));
+  run_test<3> (test3d_2, &(index[1]));
+
+  deallog << "Testing Order 3" << std::endl;
+  run_test<2> (test2d_3, &(index[2]));
+  run_test<3> (test3d_3, &(index[2]));
+      
+  return 0;
+}

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.