]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Made degree of quadrature depending on the degree of the finite element used in the...
authorRoland <roland.richter@ntnu.no>
Thu, 25 Apr 2019 15:06:20 +0000 (17:06 +0200)
committerRoland <roland.richter@ntnu.no>
Thu, 25 Apr 2019 15:06:36 +0000 (17:06 +0200)
25 files changed:
doc/news/changes/major/20190425RolandRichter [new file with mode: 0644]
examples/step-13/step-13.cc
examples/step-14/step-14.cc
examples/step-15/step-15.cc
examples/step-16b/doc/results.dox
examples/step-16b/step-16b.cc
examples/step-17/step-17.cc
examples/step-18/step-18.cc
examples/step-23/step-23.cc
examples/step-24/step-24.cc
examples/step-25/step-25.cc
examples/step-29/step-29.cc
examples/step-3/doc/results.dox
examples/step-3/step-3.cc
examples/step-33/step-33.cc
examples/step-36/step-36.cc
examples/step-4/step-4.cc
examples/step-40/step-40.cc
examples/step-5/step-5.cc
examples/step-51/step-51.cc
examples/step-55/step-55.cc
examples/step-6/step-6.cc
examples/step-7/doc/results.dox
examples/step-7/step-7.cc
examples/step-8/step-8.cc

diff --git a/doc/news/changes/major/20190425RolandRichter b/doc/news/changes/major/20190425RolandRichter
new file mode 100644 (file)
index 0000000..7477577
--- /dev/null
@@ -0,0 +1,3 @@
+New: The tutorial examples now all use an interpolation formula with a degree depending on the degree of the finite element, which can be changed dynamically.
+<br>
+(Roland Richter, 2019/04/25)
index b37aac8107ab1979c2da863305fe8d194e6a61a8..71613bf13e1f1013c6c33d876fd95198a5d52f4e 100644 (file)
@@ -1221,7 +1221,7 @@ namespace Step13
         this->triangulation->n_active_cells());
       KellyErrorEstimator<dim>::estimate(
         this->dof_handler,
-        QGauss<dim - 1>(3),
+        QGauss<dim - 1>(this->fe->degree + 1),
         std::map<types::boundary_id, const Function<dim> *>(),
         this->solution,
         estimated_error_per_cell);
index b79fef2d2f80a4f7e8643ee9f269c8cf5729d976..6e1c319a563c8dbd81fb3244958931561da5af2a 100644 (file)
@@ -878,7 +878,7 @@ namespace Step14
         this->triangulation->n_active_cells());
       KellyErrorEstimator<dim>::estimate(
         this->dof_handler,
-        QGauss<dim - 1>(3),
+        QGauss<dim - 1>(this->fe->degree + 1),
         std::map<types::boundary_id, const Function<dim> *>(),
         this->solution,
         estimated_error_per_cell);
@@ -1530,7 +1530,7 @@ namespace Step14
       // Initialize a <code>FEValues</code> object with a quadrature formula,
       // have abbreviations for the number of quadrature points and shape
       // functions...
-      QGauss<dim>        quadrature(4);
+      QGauss<dim>        quadrature(dof_handler.get_fe().degree + 1);
       FEValues<dim>      fe_values(dof_handler.get_fe(),
                               quadrature,
                               update_gradients | update_quadrature_points |
index 0e980c9e7efd23132781ed994f9e0f00179f9eb2..8f5c01b9ef90a44b0fd05ad0eb01985094a14256 100644 (file)
@@ -234,7 +234,7 @@ namespace Step15
   template <int dim>
   void MinimalSurfaceProblem<dim>::assemble_system()
   {
-    const QGauss<dim> quadrature_formula(3);
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
 
     system_matrix = 0;
     system_rhs    = 0;
@@ -376,7 +376,7 @@ namespace Step15
 
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(3),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       present_solution,
       estimated_error_per_cell);
@@ -504,7 +504,7 @@ namespace Step15
     evaluation_point = present_solution;
     evaluation_point.add(alpha, newton_update);
 
-    const QGauss<dim> quadrature_formula(3);
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
     FEValues<dim>     fe_values(fe,
                             quadrature_formula,
                             update_gradients | update_quadrature_points |
index 4f5baa31f477fece7dc26ceb0b26c3fb5619811f..cd4951e1daf23ef9390675c6bf000fde0ce11578 100644 (file)
@@ -15,38 +15,38 @@ DEAL::   Number of degrees of freedom: 25 (by level: 8, 25)
 DEAL:cg::Starting value 0.510691
 DEAL:cg::Convergence step 6 value 4.59193e-14
 DEAL::Cycle 1
-DEAL::   Number of active cells:       41
-DEAL::   Number of degrees of freedom: 52 (by level: 8, 25, 41)
-DEAL:cg::Starting value 0.455356
-DEAL:cg::Convergence step 8 value 3.09682e-13
+DEAL::   Number of active cells:       44
+DEAL::   Number of degrees of freedom: 55 (by level: 8, 25, 45)
+DEAL:cg::Starting value 0.440678
+DEAL:cg::Convergence step 8 value 1.99419e-13
 DEAL::Cycle 2
-DEAL::   Number of active cells:       80
-DEAL::   Number of degrees of freedom: 100 (by level: 8, 25, 61, 52)
-DEAL:cg::Starting value 0.394469
-DEAL:cg::Convergence step 9 value 1.96993e-13
+DEAL::   Number of active cells:       86
+DEAL::   Number of degrees of freedom: 105 (by level: 8, 25, 69, 49)
+DEAL:cg::Starting value 0.371855
+DEAL:cg::Convergence step 9 value 1.13984e-13
 DEAL::Cycle 3
-DEAL::   Number of active cells:       161
-DEAL::   Number of degrees of freedom: 190 (by level: 8, 25, 77, 160)
-DEAL:cg::Starting value 0.322156
-DEAL:cg::Convergence step 9 value 2.94418e-13
+DEAL::   Number of active cells:       170
+DEAL::   Number of degrees of freedom: 200 (by level: 8, 25, 77, 174)
+DEAL:cg::Starting value 0.318967
+DEAL:cg::Convergence step 9 value 2.62112e-13
 DEAL::Cycle 4
-DEAL::   Number of active cells:       311
-DEAL::   Number of degrees of freedom: 364 (by level: 8, 25, 86, 227, 174)
-DEAL:cg::Starting value 0.279667
-DEAL:cg::Convergence step 10 value 3.45746e-13
+DEAL::   Number of active cells:       332
+DEAL::   Number of degrees of freedom: 388 (by level: 8, 25, 86, 231, 204)
+DEAL:cg::Starting value 0.276534
+DEAL:cg::Convergence step 10 value 1.69562e-13
 DEAL::Cycle 5
-DEAL::   Number of active cells:       593
-DEAL::   Number of degrees of freedom: 667 (by level: 8, 25, 89, 231, 490, 96)
-DEAL:cg::Starting value 0.215917
-DEAL:cg::Convergence step 10 value 1.03758e-13
+DEAL::   Number of active cells:       632
+DEAL::   Number of degrees of freedom: 714 (by level: 8, 25, 89, 231, 514, 141)
+DEAL:cg::Starting value 0.215300
+DEAL:cg::Convergence step 10 value 6.47463e-13
 DEAL::Cycle 6
-DEAL::   Number of active cells:       1127
-DEAL::   Number of degrees of freedom: 1251 (by level: 8, 25, 89, 274, 760, 417, 178)
-DEAL:cg::Starting value 0.185906
-DEAL:cg::Convergence step 10 value 3.40351e-13
+DEAL::   Number of active cells:       1202
+DEAL::   Number of degrees of freedom: 1332 (by level: 8, 25, 89, 282, 771, 435, 257)
+DEAL:cg::Starting value 0.175848
+DEAL:cg::Convergence step 10 value 1.80664e-13
 DEAL::Cycle 7
-DEAL::   Number of active cells:       2144
-DEAL::   Number of degrees of freedom: 2359 (by level: 8, 25, 89, 308, 779, 1262, 817)
-DEAL:cg::Starting value 0.141519
-DEAL:cg::Convergence step 10 value 5.74965e-13
+DEAL::   Number of active cells:       2288
+DEAL::   Number of degrees of freedom: 2511 (by level: 8, 25, 89, 318, 779, 1420, 829, 30)
+DEAL:cg::Starting value 0.136724
+DEAL:cg::Convergence step 11 value 9.73331e-14
 </pre>
index b35c6d4f558b620bbbe2f0c02f6e9ba8046d6023..101a570c4678a9130114f643d2f0d0810c897a51 100644 (file)
@@ -561,7 +561,7 @@ namespace Step16
 
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(3),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       solution,
       estimated_error_per_cell);
index c837d09506bea99c998d4b2361377b49a1e1a91d..14cd33023bd2f12c1becb47814bd132979d4bfb6 100644 (file)
@@ -462,7 +462,7 @@ namespace Step17
   template <int dim>
   void ElasticProblem<dim>::assemble_system()
   {
-    QGauss<dim>   quadrature_formula(2);
+    QGauss<dim>   quadrature_formula(fe.degree + 1);
     FEValues<dim> fe_values(fe,
                             quadrature_formula,
                             update_values | update_gradients |
@@ -765,7 +765,7 @@ namespace Step17
     Vector<float> local_error_per_cell(triangulation.n_active_cells());
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(2),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       localized_solution,
       local_error_per_cell,
index 5796801b572fe8a59baa517e3878ae1e077d4a84..17a231225aaf7d7b603c992be5ee93964da9b830 100644 (file)
@@ -716,7 +716,7 @@ namespace Step18
     : triangulation(MPI_COMM_WORLD)
     , fe(FE_Q<dim>(1), dim)
     , dof_handler(triangulation)
-    , quadrature_formula(2)
+    , quadrature_formula(fe.degree + 1)
     , present_time(0.0)
     , present_timestep(1.0)
     , end_time(10.0)
@@ -1441,7 +1441,7 @@ namespace Step18
     Vector<float> error_per_cell(triangulation.n_active_cells());
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(2),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       incremental_displacement,
       error_per_cell,
index a120d0991fd68b652fa4ba62bb756d8877473e21..5b814445501b6377e9c59743c231a55ab4fcaf9e 100644 (file)
@@ -373,9 +373,11 @@ namespace Step23
     matrix_u.reinit(sparsity_pattern);
     matrix_v.reinit(sparsity_pattern);
 
-    MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+    MatrixCreator::create_mass_matrix(dof_handler,
+                                      QGauss<dim>(fe.degree + 1),
+                                      mass_matrix);
     MatrixCreator::create_laplace_matrix(dof_handler,
-                                         QGauss<dim>(3),
+                                         QGauss<dim>(fe.degree + 1),
                                          laplace_matrix);
 
     // The rest of the function is spent on setting vector sizes to the
@@ -478,12 +480,12 @@ namespace Step23
 
     VectorTools::project(dof_handler,
                          constraints,
-                         QGauss<dim>(3),
+                         QGauss<dim>(fe.degree + 1),
                          InitialValuesU<dim>(),
                          old_solution_u);
     VectorTools::project(dof_handler,
                          constraints,
-                         QGauss<dim>(3),
+                         QGauss<dim>(fe.degree + 1),
                          InitialValuesV<dim>(),
                          old_solution_v);
 
index 737d5cd5537627cfdb9b4bfe48b1f47266a28f3f..ed0622a618ad28b6724ffffa3e3fbda0d2693020 100644 (file)
@@ -289,9 +289,11 @@ namespace Step24
     mass_matrix.reinit(sparsity_pattern);
     laplace_matrix.reinit(sparsity_pattern);
 
-    MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+    MatrixCreator::create_mass_matrix(dof_handler,
+                                      QGauss<dim>(fe.degree + 1),
+                                      mass_matrix);
     MatrixCreator::create_laplace_matrix(dof_handler,
-                                         QGauss<dim>(3),
+                                         QGauss<dim>(fe.degree + 1),
                                          laplace_matrix);
 
     // The second difference, as mentioned, to step-23 is that we need to
@@ -336,7 +338,7 @@ namespace Step24
     // do something only if that particular face is at the boundary of the
     // domain. Like this:
     {
-      const QGauss<dim - 1> quadrature_formula(3);
+      const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
       FEFaceValues<dim>     fe_values(fe,
                                   quadrature_formula,
                                   update_values | update_JxW_values);
@@ -472,7 +474,7 @@ namespace Step24
 
     VectorTools::project(dof_handler,
                          constraints,
-                         QGauss<dim>(3),
+                         QGauss<dim>(fe.degree + 1),
                          InitialValuesP<dim>(),
                          old_solution_p);
     old_solution_v = 0;
index 8f8d4d7b21bc660c6a02caec3b289daa903ad6eb..7372a4798003a8d8f5699bbbfd9c0ffab01f8a53 100644 (file)
@@ -308,9 +308,11 @@ namespace Step25
     mass_matrix.reinit(sparsity_pattern);
     laplace_matrix.reinit(sparsity_pattern);
 
-    MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+    MatrixCreator::create_mass_matrix(dof_handler,
+                                      QGauss<dim>(fe.degree + 1),
+                                      mass_matrix);
     MatrixCreator::create_laplace_matrix(dof_handler,
-                                         QGauss<dim>(3),
+                                         QGauss<dim>(fe.degree + 1),
                                          laplace_matrix);
 
     solution.reinit(dof_handler.n_dofs());
@@ -398,7 +400,7 @@ namespace Step25
                                                Vector<double> &nl_term) const
   {
     nl_term = 0;
-    const QGauss<dim> quadrature_formula(3);
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
     FEValues<dim>     fe_values(fe,
                             quadrature_formula,
                             update_values | update_JxW_values |
@@ -462,7 +464,7 @@ namespace Step25
     const Vector<double> &new_data,
     SparseMatrix<double> &nl_matrix) const
   {
-    QGauss<dim>   quadrature_formula(3);
+    QGauss<dim>   quadrature_formula(fe.degree + 1);
     FEValues<dim> fe_values(fe,
                             quadrature_formula,
                             update_values | update_JxW_values |
@@ -601,7 +603,7 @@ namespace Step25
       constraints.close();
       VectorTools::project(dof_handler,
                            constraints,
-                           QGauss<dim>(3),
+                           QGauss<dim>(fe.degree + 1),
                            InitialValues<dim>(1, time),
                            solution);
     }
index cc2d2b4590376fdf15037ca2c9f39b724563035f..805f55801f9576bdee20c9fe7d4456270a439e86 100644 (file)
@@ -554,8 +554,8 @@ namespace Step29
     // used. Since our bilinear form involves boundary integrals on
     // $\Gamma_2$, we also need a quadrature rule for surface integration on
     // the faces, which are $dim-1$ dimensional:
-    QGauss<dim>     quadrature_formula(2);
-    QGauss<dim - 1> face_quadrature_formula(2);
+    QGauss<dim>     quadrature_formula(fe.degree + 1);
+    QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 
     const unsigned int n_q_points      = quadrature_formula.size(),
                        n_face_q_points = face_quadrature_formula.size(),
index a39a00addc961d2a7f4af374b2d55ff1b6a91ff5..6e61813cf71aaa9642fae180c9b2195c92e6bce4 100644 (file)
@@ -192,7 +192,7 @@ suggestions:
   @code
     std::cout << "Mean value: "
               << VectorTools::compute_mean_value (dof_handler,
-                                                 QGauss<2>(3),
+                                                 QGauss<2>(fe.degree + 1),
                                                  solution,
                                                  0)
               << std::endl;
index 1e50eef4db8dfef605654429d6f27bf9b5840808..9f35b732793c61694fe44c0a5545df9803a14d98 100644 (file)
@@ -272,7 +272,7 @@ void Step3::assemble_system()
   // 2D. This quadrature formula integrates polynomials of degrees up to three
   // exactly (in 1D). It is easy to check that this is sufficient for the
   // present problem:
-  QGauss<2> quadrature_formula(2);
+  QGauss<2> quadrature_formula(fe.degree + 1);
   // And we initialize the object which we have briefly talked about above. It
   // needs to be told which finite element we want to use, and the quadrature
   // points and their weights (jointly described by a Quadrature object). As
index 917b51fcd30db0b143390d4cc6adc5f75a2b83ca..a24a5665fc25e78cdbd04858d640bcd4b680edbb 100644 (file)
@@ -1381,8 +1381,8 @@ namespace Step33
     : mapping()
     , fe(FE_Q<dim>(1), EulerEquations<dim>::n_components)
     , dof_handler(triangulation)
-    , quadrature(2)
-    , face_quadrature(2)
+    , quadrature(fe.degree + 1)
+    , face_quadrature(fe.degree + 1)
     , verbose_cout(std::cout, false)
   {
     ParameterHandler prm;
index fe6e593e15d10319fd0701b52ebf1e32e1af475d..6213053568c4cb1739d2592d7450ade7b390f678 100644 (file)
@@ -237,7 +237,7 @@ namespace Step36
   template <int dim>
   void EigenvalueProblem<dim>::assemble_system()
   {
-    QGauss<dim> quadrature_formula(2);
+    QGauss<dim> quadrature_formula(fe.degree + 1);
 
     FEValues<dim> fe_values(fe,
                             quadrature_formula,
index a4b0aa53907b8fc30b21ef44bcd26f060ce8afc9..e3d1447b228865f359081aa7cf35609e161bf36c 100644 (file)
@@ -308,7 +308,7 @@ void Step4<dim>::setup_system()
 template <int dim>
 void Step4<dim>::assemble_system()
 {
-  QGauss<dim> quadrature_formula(2);
+  QGauss<dim> quadrature_formula(fe.degree + 1);
 
   // We wanted to have a non-constant right hand side, so we use an object of
   // the class declared above to generate the necessary data. Since this right
index 0837dd65f6623169aad427920de31e3e1609575a..25837b4726bae817742e5f282668943ecdea7072 100644 (file)
@@ -366,7 +366,7 @@ namespace Step40
   {
     TimerOutput::Scope t(computing_timer, "assembly");
 
-    const QGauss<dim> quadrature_formula(3);
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
 
     FEValues<dim> fe_values(fe,
                             quadrature_formula,
@@ -523,7 +523,7 @@ namespace Step40
     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(3),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       locally_relevant_solution,
       estimated_error_per_cell);
index 8b694d173a4f2e65acb09f8fc6cfbbdd63186ba9..2c0af87460693ec758e6ab9b8ae44865587c254b 100644 (file)
@@ -162,7 +162,7 @@ void Step5<dim>::setup_system()
 template <int dim>
 void Step5<dim>::assemble_system()
 {
-  QGauss<dim> quadrature_formula(2);
+  QGauss<dim> quadrature_formula(fe.degree + 1);
 
   FEValues<dim> fe_values(fe,
                           quadrature_formula,
index c44dbcda51da41de5137a84a417e598b976ca48e..36a1d91145fa0c6cf969153191098fb582bf02e4 100644 (file)
@@ -1333,7 +1333,7 @@ namespace Step51
               std::map<types::boundary_id, const Function<dim> *>
                 neumann_boundary;
               KellyErrorEstimator<dim>::estimate(dof_handler_local,
-                                                 QGauss<dim - 1>(3),
+                                                 QGauss<dim - 1>(fe.degree + 1),
                                                  neumann_boundary,
                                                  solution_local,
                                                  estimated_error_per_cell,
index 2d42ef85f43234790cd6550d94527f88e69b7fbc..10b6ac513f102b1db93fd6b76a37fddfabd10910 100644 (file)
@@ -691,7 +691,7 @@ namespace Step55
         FEValuesExtractors::Vector velocities(0);
         KellyErrorEstimator<dim>::estimate(
           dof_handler,
-          QGauss<dim - 1>(3),
+          QGauss<dim - 1>(fe.degree + 1),
           std::map<types::boundary_id, const Function<dim> *>(),
           locally_relevant_solution,
           estimated_error_per_cell,
index b45caaf7bd4f15af8104c460f49c5401d6af5e99..18bb3b476823c404fcc17dbb6d04c694c3bcbcbd 100644 (file)
@@ -269,7 +269,7 @@ void Step6<dim>::setup_system()
 template <int dim>
 void Step6<dim>::assemble_system()
 {
-  const QGauss<dim> quadrature_formula(3);
+  const QGauss<dim> quadrature_formula(fe.degree + 1);
 
   FEValues<dim> fe_values(fe,
                           quadrature_formula,
@@ -422,7 +422,7 @@ void Step6<dim>::refine_grid()
 
   KellyErrorEstimator<dim>::estimate(
     dof_handler,
-    QGauss<dim - 1>(3),
+    QGauss<dim - 1>(fe.degree + 1),
     std::map<types::boundary_id, const Function<dim> *>(),
     solution,
     estimated_error_per_cell);
index 6671886260094b2a84041113ba896beeeeebe13f..fdc5e9a67b084a91b9ad634220493d42d70058ff 100644 (file)
@@ -53,15 +53,15 @@ Cycle 8:
    Number of degrees of freedom: 15912
 
 cycle cells dofs     L2        H1      Linfty   
-    0    64    81 1.576e-01 1.418e+00 2.707e-01 
-    1   124   157 4.285e-02 1.285e+00 1.469e-01 
-    2   280   341 1.593e-02 7.909e-01 8.034e-02 
-    3   577   690 9.359e-03 5.096e-01 2.784e-02 
-    4  1099  1264 2.865e-03 3.038e-01 9.822e-03 
-    5  2191  2452 1.480e-03 2.106e-01 5.679e-03 
-    6  4165  4510 6.907e-04 1.462e-01 2.338e-03 
-    7  7915  8440 4.743e-04 1.055e-01 1.442e-03 
-    8 15196 15912 1.920e-04 7.468e-02 7.259e-04 
+    0    64    81 1.840e+00 2.858e+00 1.835e+00 
+    1   124   157 5.190e-02 1.200e+00 1.344e-01 
+    2   280   341 1.439e-02 7.892e-01 7.554e-02 
+    3   577   690 8.627e-03 5.061e-01 2.805e-02 
+    4  1099  1264 3.217e-03 3.030e-01 1.073e-02 
+    5  2191  2452 1.445e-03 2.097e-01 5.073e-03 
+    6  4165  4510 8.387e-04 1.460e-01 2.013e-03 
+    7  7915  8440 7.051e-04 1.053e-01 1.804e-03 
+    8 15196 15912 2.774e-04 7.463e-02 6.911e-04 
 
 Solving with Q1 elements, global refinement
 ===========================================
@@ -83,18 +83,18 @@ Cycle 4:
    Number of degrees of freedom: 16641
 
 cycle cells dofs     L2        H1      Linfty   
-    0    64    81 1.576e-01 1.418e+00 2.707e-01 
-    1   256   289 4.280e-02 1.285e+00 1.444e-01 
-    2  1024  1089 1.352e-02 7.556e-01 7.772e-02 
-    3  4096  4225 3.423e-03 3.822e-01 2.332e-02 
-    4 16384 16641 8.586e-04 1.917e-01 6.097e-03 
-
-n cells         H1                  L2          
-0    64 1.418e+00    -    - 1.576e-01    -    - 
-1   256 1.285e+00 1.10 0.14 4.280e-02 3.68 1.88 
-2  1024 7.556e-01 1.70 0.77 1.352e-02 3.17 1.66 
-3  4096 3.822e-01 1.98 0.98 3.423e-03 3.95 1.98 
-4 16384 1.917e-01 1.99 1.00 8.586e-04 3.99 2.00 
+    0    64    81 1.840e+00 2.858e+00 1.835e+00 
+    1   256   289 3.570e-02 1.199e+00 1.307e-01 
+    2  1024  1089 1.192e-02 7.565e-01 7.168e-02 
+    3  4096  4225 3.047e-03 3.823e-01 2.128e-02 
+    4 16384 16641 7.660e-04 1.917e-01 5.554e-03 
+
+n cells         H1                   L2          
+0    64 2.858e+00    -    - 1.840e+00     -    - 
+1   256 1.199e+00 2.38 1.25 3.570e-02 51.54 5.69 
+2  1024 7.565e-01 1.58 0.66 1.192e-02  2.99 1.58 
+3  4096 3.823e-01 1.98 0.98 3.047e-03  3.91 1.97 
+4 16384 1.917e-01 1.99 1.00 7.660e-04  3.98 1.99 
 
 Solving with Q2 elements, global refinement
 ===========================================
@@ -169,7 +169,7 @@ cycle cells dofs     L2        H1      Linfty
     5  2059  9223 7.738e-05 1.974e-02 7.270e-04 
     6  3913 17887 2.925e-05 8.772e-03 1.463e-04 
     7  7441 33807 1.024e-05 4.121e-03 8.567e-05 
-    8 14212 64731 3.761e-06 2.108e-03 2.167e-05 
+    8 14212 64731 3.761e-06 2.108e-03 2.167e-05  
 @endcode
 
 
@@ -192,11 +192,8 @@ here).
 
 Go ahead and run the program with higher order elements (Q3, Q4, ...). You
 will notice that assertions in several parts of the code will trigger (for
-example in the generation of the filename for the data output). After fixing
-these you will not see the correct convergence orders that the theory
-predicts. This is because the orders for the quadrature formulas are
-hard-coded in this program and this order is not enough for higher order
-discretizations. What is a good way to pick the orders dynamically?
+example in the generation of the filename for the data output). You might have to address these,
+but it should not be very hard to get the program to work!
 
 <h4> Convergence Comparison </h4>
 
index 97c0e7daa91c9180b20fe06dace41b38fa6f67ee..223460cf7dcfe20bf95230d82dc7780b44a75d09 100644 (file)
@@ -541,8 +541,8 @@ namespace Step7
   template <int dim>
   void HelmholtzProblem<dim>::assemble_system()
   {
-    QGauss<dim>     quadrature_formula(3);
-    QGauss<dim - 1> face_quadrature_formula(3);
+    QGauss<dim>     quadrature_formula(fe->degree + 1);
+    QGauss<dim - 1> face_quadrature_formula(fe->degree + 1);
 
     const unsigned int n_q_points      = quadrature_formula.size();
     const unsigned int n_face_q_points = face_quadrature_formula.size();
@@ -797,7 +797,7 @@ namespace Step7
 
             KellyErrorEstimator<dim>::estimate(
               dof_handler,
-              QGauss<dim - 1>(3),
+              QGauss<dim - 1>(fe->degree + 1),
               std::map<types::boundary_id, const Function<dim> *>(),
               solution,
               estimated_error_per_cell);
@@ -853,7 +853,7 @@ namespace Step7
                                       solution,
                                       Solution<dim>(),
                                       difference_per_cell,
-                                      QGauss<dim>(3),
+                                      QGauss<dim>(fe->degree + 1),
                                       VectorTools::L2_norm);
     const double L2_error =
       VectorTools::compute_global_error(triangulation,
@@ -871,7 +871,7 @@ namespace Step7
                                       solution,
                                       Solution<dim>(),
                                       difference_per_cell,
-                                      QGauss<dim>(3),
+                                      QGauss<dim>(fe->degree + 1),
                                       VectorTools::H1_seminorm);
     const double H1_error =
       VectorTools::compute_global_error(triangulation,
@@ -883,8 +883,9 @@ namespace Step7
     // points. Since this depends quite sensitively on the quadrature rule
     // being used, and since we would like to avoid false results due to
     // super-convergence effects at some points, we use a special quadrature
-    // rule that is obtained by iterating the trapezoidal rule five times in
-    // each space direction. Note that the constructor of the QIterated class
+    // rule that is obtained by iterating the trapezoidal rule by the degree of
+    // of the finite element times two plus one in each space direction.
+    // Note that the constructor of the QIterated class
     // takes a one-dimensional quadrature rule and a number that tells it how
     // often it shall use this rule in each space direction.
     //
@@ -893,7 +894,7 @@ namespace Step7
     // from the L infinity errors on each cell with a call to
     // VectorTools::compute_global_error.
     const QTrapez<1>     q_trapez;
-    const QIterated<dim> q_iterated(q_trapez, 5);
+    const QIterated<dim> q_iterated(q_trapez, fe->degree * 2 + 1);
     VectorTools::integrate_difference(dof_handler,
                                       solution,
                                       Solution<dim>(),
index 4927b72858b46d6cbe44790c9c2f9d13be2902e0..3f7318fd21c735c410ea48380c35bcd0469a89a4 100644 (file)
@@ -269,7 +269,7 @@ namespace Step8
   template <int dim>
   void ElasticProblem<dim>::assemble_system()
   {
-    QGauss<dim> quadrature_formula(2);
+    QGauss<dim> quadrature_formula(fe.degree + 1);
 
     FEValues<dim> fe_values(fe,
                             quadrature_formula,
@@ -490,7 +490,7 @@ namespace Step8
 
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(2),
+      QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(),
       solution,
       estimated_error_per_cell);

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.