Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-74.h
Go to the documentation of this file.
1
274) const
275 *   {
276 *   using numbers::PI;
277 *   for (unsigned int i = 0; i < values.size(); ++i)
278 *   values[i] =
279 *   std::sin(2. * PI * points[i][0]) * std::sin(2. * PI * points[i][1]);
280 *   }
281 *  
282 *  
283 *  
284 *   template <int dim>
285 *   Tensor<1, dim>
286 *   SmoothSolution<dim>::gradient(const Point<dim> &point,
287 *   const unsigned int /*component*/) const
288 *   {
289 *   Tensor<1, dim> return_value;
290 *   using numbers::PI;
291 *   return_value[0] =
292 *   2. * PI * std::cos(2. * PI * point[0]) * std::sin(2. * PI * point[1]);
293 *   return_value[1] =
294 *   2. * PI * std::sin(2. * PI * point[0]) * std::cos(2. * PI * point[1]);
295 *   return return_value;
296 *   }
297 *  
298 *  
299 *  
300 * @endcode
301 *
302 * The corresponding right-hand side of the smooth function:
303 *
304 * @code
305 *   template <int dim>
306 *   class SmoothRightHandSide : public Function<dim>
307 *   {
308 *   public:
309 *   SmoothRightHandSide()
310 *   : Function<dim>()
311 *   {}
312 *  
313 *   virtual void value_list(const std::vector<Point<dim>> &points,
314 *   std::vector<double> &values,
315 *   const unsigned int /*component*/) const override;
316 *   };
317 *  
318 *  
319 *  
320 *   template <int dim>
321 *   void
322 *   SmoothRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
323 *   std::vector<double> &values,
324 *   const unsigned int /*component*/) const
325 *   {
326 *   using numbers::PI;
327 *   for (unsigned int i = 0; i < values.size(); ++i)
328 *   values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) *
329 *   std::sin(2. * PI * points[i][1]);
330 *   }
331 *  
332 *  
333 *  
334 * @endcode
335 *
336 * The right-hand side that corresponds to the function
338 * assume that the diffusion coefficient @f$\nu = 1@f$:
339 *
340 * @code
341 *   template <int dim>
342 *   class SingularRightHandSide : public Function<dim>
343 *   {
344 *   public:
345 *   SingularRightHandSide()
346 *   : Function<dim>()
347 *   {}
348 *  
349 *   virtual void value_list(const std::vector<Point<dim>> &points,
350 *   std::vector<double> &values,
351 *   const unsigned int /*component*/) const override;
352 *  
353 *   private:
354 *   const Functions::LSingularityFunction ref;
355 *   };
356 *  
357 *  
358 *  
359 *   template <int dim>
360 *   void
361 *   SingularRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
362 *   std::vector<double> &values,
363 *   const unsigned int /*component*/) const
364 *   {
365 *   for (unsigned int i = 0; i < values.size(); ++i)
366 *   values[i] = -ref.laplacian(points[i]);
367 *   }
368 *  
369 *  
370 *  
371 * @endcode
372 *
373 *
374 * <a name="step_74-Auxiliaryfunctions"></a>
375 * <h3>Auxiliary functions</h3>
376 * This function computes the penalty @f$\sigma@f$.
377 *
378 * @code
379 *   double get_penalty_factor(const unsigned int fe_degree,
380 *   const double cell_extent_left,
381 *   const double cell_extent_right)
382 *   {
383 *   const unsigned int degree = std::max(1U, fe_degree);
384 *   return degree * (degree + 1.) * 0.5 *
385 *   (1. / cell_extent_left + 1. / cell_extent_right);
386 *   }
387 *  
388 *  
389 * @endcode
390 *
391 *
392 * <a name="step_74-TheCopyData"></a>
393 * <h3>The CopyData</h3>
394 * In the following, we define "Copy" objects for the MeshWorker::mesh_loop(),
395 * which is essentially the same as @ref step_12 "step-12". Note that the
396 * "Scratch" object is not defined here because we use
397 * MeshWorker::ScratchData<dim> instead. (The use of "Copy" and "Scratch"
398 * objects is extensively explained in the WorkStream namespace documentation.
399 *
400 * @code
401 *   struct CopyDataFace
402 *   {
404 *   std::vector<types::global_dof_index> joint_dof_indices;
405 *   std::array<double, 2> values;
406 *   std::array<unsigned int, 2> cell_indices;
407 *   };
408 *  
409 *  
410 *  
411 *   struct CopyData
412 *   {
414 *   Vector<double> cell_rhs;
415 *   std::vector<types::global_dof_index> local_dof_indices;
416 *   std::vector<CopyDataFace> face_data;
417 *   double value;
418 *   unsigned int cell_index;
419 *  
420 *  
421 *   template <class Iterator>
422 *   void reinit(const Iterator &cell, const unsigned int dofs_per_cell)
423 *   {
424 *   cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
425 *   cell_rhs.reinit(dofs_per_cell);
426 *   local_dof_indices.resize(dofs_per_cell);
427 *   cell->get_dof_indices(local_dof_indices);
428 *   }
429 *   };
430 *  
431 *  
432 *  
433 * @endcode
434 *
435 *
436 * <a name="step_74-TheSIPGLaplaceclass"></a>
437 * <h3>The SIPGLaplace class</h3>
438 * After these preparations, we proceed with the main class of this program,
439 * called `SIPGLaplace`. The overall structure of the class is as in many
440 * of the other tutorial programs. Major differences will only come up in the
441 * implementation of the assemble functions, since we use FEInterfaceValues to
442 * assemble face terms.
443 *
444 * @code
445 *   template <int dim>
446 *   class SIPGLaplace
447 *   {
448 *   public:
449 *   SIPGLaplace(const TestCase &test_case);
450 *   void run();
451 *  
452 *   private:
453 *   void setup_system();
454 *   void assemble_system();
455 *   void solve();
456 *   void refine_grid();
457 *   void output_results(const unsigned int cycle) const;
458 *  
459 *   void compute_errors();
460 *   void compute_error_estimate();
461 *   double compute_energy_norm_error();
462 *  
464 *   const unsigned int degree;
465 *   const QGauss<dim> quadrature;
466 *   const QGauss<dim - 1> face_quadrature;
467 *   const QGauss<dim> quadrature_overintegration;
468 *   const QGauss<dim - 1> face_quadrature_overintegration;
469 *   const MappingQ1<dim> mapping;
470 *  
471 *   using ScratchData = MeshWorker::ScratchData<dim>;
472 *  
473 *   const FE_DGQ<dim> fe;
474 *   DoFHandler<dim> dof_handler;
475 *  
476 *   SparsityPattern sparsity_pattern;
477 *   SparseMatrix<double> system_matrix;
478 *   Vector<double> solution;
479 *   Vector<double> system_rhs;
480 *  
481 * @endcode
482 *
483 * The remainder of the class's members are used for the following:
484 * - Vectors to store error estimator square and energy norm square per
485 * cell.
486 * - Print convergence rate and errors on the screen.
487 * - The fiffusion coefficient @f$\nu@f$ is set to 1.
488 * - Members that store information about the test case to be computed.
489 *
490 * @code
491 *   Vector<double> estimated_error_square_per_cell;
492 *   Vector<double> energy_norm_square_per_cell;
493 *  
494 *   ConvergenceTable convergence_table;
495 *  
496 *   const double diffusion_coefficient = 1.;
497 *  
498 *   const TestCase test_case;
499 *   std::unique_ptr<const Function<dim>> exact_solution;
500 *   std::unique_ptr<const Function<dim>> rhs_function;
501 *   };
502 *  
503 * @endcode
504 *
505 * The constructor here takes the test case as input and then
506 * determines the correct solution and right-hand side classes. The
507 * remaining member variables are initialized in the obvious way.
508 *
509 * @code
510 *   template <int dim>
511 *   SIPGLaplace<dim>::SIPGLaplace(const TestCase &test_case)
512 *   : degree(3)
513 *   , quadrature(degree + 1)
514 *   , face_quadrature(degree + 1)
515 *   , quadrature_overintegration(degree + 2)
516 *   , face_quadrature_overintegration(degree + 2)
517 *   , mapping()
518 *   , fe(degree)
519 *   , dof_handler(triangulation)
520 *   , test_case(test_case)
521 *   {
522 *   if (test_case == TestCase::convergence_rate)
523 *   {
524 *   exact_solution = std::make_unique<const SmoothSolution<dim>>();
525 *   rhs_function = std::make_unique<const SmoothRightHandSide<dim>>();
526 *   }
527 *  
528 *   else if (test_case == TestCase::l_singularity)
529 *   {
530 *   exact_solution =
531 *   std::make_unique<const Functions::LSingularityFunction>();
532 *   rhs_function = std::make_unique<const SingularRightHandSide<dim>>();
533 *   }
534 *   else
535 *   AssertThrow(false, ExcNotImplemented());
536 *   }
537 *  
538 *  
539 *  
540 *   template <int dim>
541 *   void SIPGLaplace<dim>::setup_system()
542 *   {
543 *   dof_handler.distribute_dofs(fe);
544 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
545 *   DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
546 *   sparsity_pattern.copy_from(dsp);
547 *  
548 *   system_matrix.reinit(sparsity_pattern);
549 *   solution.reinit(dof_handler.n_dofs());
550 *   system_rhs.reinit(dof_handler.n_dofs());
551 *   }
552 *  
553 *  
554 *  
555 * @endcode
556 *
557 *
558 * <a name="step_74-Theassemble_systemfunction"></a>
559 * <h3>The assemble_system function</h3>
560 * The assemble function here is similar to that in @ref step_12 "step-12" and @ref step_47 "step-47".
561 * Different from assembling by hand, we just need to focus
562 * on assembling on each cell, each boundary face, and each
563 * interior face. The loops over cells and faces are handled
564 * automatically by MeshWorker::mesh_loop().
565 *
566
567 *
568 * The function starts by defining a local (lambda) function that is
569 * used to integrate the cell terms:
570 *
571 * @code
572 *   template <int dim>
573 *   void SIPGLaplace<dim>::assemble_system()
574 *   {
575 *   const auto cell_worker =
576 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
577 *   ScratchData &scratch_data,
578 *   CopyData &copy_data) {
579 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
580 *   const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
581 *   copy_data.reinit(cell, dofs_per_cell);
582 *  
583 *   const std::vector<Point<dim>> &q_points =
584 *   scratch_data.get_quadrature_points();
585 *   const unsigned int n_q_points = q_points.size();
586 *   const std::vector<double> &JxW = scratch_data.get_JxW_values();
587 *  
588 *   std::vector<double> rhs(n_q_points);
589 *   rhs_function->value_list(q_points, rhs);
590 *  
591 *   for (unsigned int point = 0; point < n_q_points; ++point)
592 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
593 *   {
594 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
595 *   copy_data.cell_matrix(i, j) +=
596 *   diffusion_coefficient * // nu
597 *   fe_v.shape_grad(i, point) * // grad v_h
598 *   fe_v.shape_grad(j, point) * // grad u_h
599 *   JxW[point]; // dx
600 *  
601 *   copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
602 *   rhs[point] * // f
603 *   JxW[point]; // dx
604 *   }
605 *   };
606 *  
607 * @endcode
608 *
609 * Next, we need a function that assembles face integrals on the boundary:
610 *
611 * @code
612 *   const auto boundary_worker =
613 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
614 *   const unsigned int &face_no,
615 *   ScratchData &scratch_data,
616 *   CopyData &copy_data) {
617 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
618 *  
619 *   const std::vector<Point<dim>> &q_points =
620 *   scratch_data.get_quadrature_points();
621 *   const unsigned int n_q_points = q_points.size();
622 *   const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
623 *  
624 *   const std::vector<double> &JxW = scratch_data.get_JxW_values();
625 *   const std::vector<Tensor<1, dim>> &normals =
626 *   scratch_data.get_normal_vectors();
627 *  
628 *   std::vector<double> g(n_q_points);
629 *   exact_solution->value_list(q_points, g);
630 *  
631 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
632 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
633 *  
634 *   for (unsigned int point = 0; point < n_q_points; ++point)
635 *   {
636 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
637 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
638 *   copy_data.cell_matrix(i, j) +=
639 *   (-diffusion_coefficient * // - nu
640 *   fe_fv.shape_value(i, point) * // v_h
641 *   (fe_fv.shape_grad(j, point) * // (grad u_h .
642 *   normals[point]) // n)
643 *  
644 *   - diffusion_coefficient * // - nu
645 *   (fe_fv.shape_grad(i, point) * // (grad v_h .
646 *   normals[point]) * // n)
647 *   fe_fv.shape_value(j, point) // u_h
648 *  
649 *   + diffusion_coefficient * penalty * // + nu sigma
650 *   fe_fv.shape_value(i, point) * // v_h
651 *   fe_fv.shape_value(j, point) // u_h
652 *  
653 *   ) *
654 *   JxW[point]; // dx
655 *  
656 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
657 *   copy_data.cell_rhs(i) +=
658 *   (-diffusion_coefficient * // - nu
659 *   (fe_fv.shape_grad(i, point) * // (grad v_h .
660 *   normals[point]) * // n)
661 *   g[point] // g
662 *  
663 *  
664 *   + diffusion_coefficient * penalty * // + nu sigma
665 *   fe_fv.shape_value(i, point) * g[point] // v_h g
666 *  
667 *   ) *
668 *   JxW[point]; // dx
669 *   }
670 *   };
671 *  
672 * @endcode
673 *
674 * Finally, a function that assembles face integrals on interior
675 * faces. To reinitialize FEInterfaceValues, we need to pass
676 * cells, face and subface indices (for adaptive refinement) to
677 * the reinit() function of FEInterfaceValues:
678 *
679 * @code
680 *   const auto face_worker =
681 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
682 *   const unsigned int &f,
683 *   const unsigned int &sf,
684 *   const typename DoFHandler<dim>::cell_iterator &ncell,
685 *   const unsigned int &nf,
686 *   const unsigned int &nsf,
687 *   ScratchData &scratch_data,
688 *   CopyData &copy_data) {
689 *   const FEInterfaceValues<dim> &fe_iv =
690 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
691 *  
692 *   copy_data.face_data.emplace_back();
693 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
694 *   const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
695 *   copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
696 *   copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
697 *  
698 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
699 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
700 *  
701 *   const double extent1 = cell->measure() / cell->face(f)->measure();
702 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
703 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
704 *  
705 *   for (const unsigned int point : fe_iv.quadrature_point_indices())
706 *   {
707 *   for (const unsigned int i : fe_iv.dof_indices())
708 *   for (const unsigned int j : fe_iv.dof_indices())
709 *   copy_data_face.cell_matrix(i, j) +=
710 *   (-diffusion_coefficient * // - nu
711 *   fe_iv.jump_in_shape_values(i, point) * // [v_h]
712 *   (fe_iv.average_of_shape_gradients(j,
713 *   point) * // ({grad u_h} .
714 *   normals[point]) // n)
715 *  
716 *   - diffusion_coefficient * // - nu
717 *   (fe_iv.average_of_shape_gradients(i,
718 *   point) * // (grad v_h .
719 *   normals[point]) * // n)
720 *   fe_iv.jump_in_shape_values(j, point) // [u_h]
721 *  
722 *   + diffusion_coefficient * penalty * // + nu sigma
723 *   fe_iv.jump_in_shape_values(i, point) * // [v_h]
724 *   fe_iv.jump_in_shape_values(j, point) // [u_h]
725 *  
726 *   ) *
727 *   JxW[point]; // dx
728 *   }
729 *   };
730 *  
731 * @endcode
732 *
733 * The following lambda function will then copy data into the
734 * global matrix and right-hand side. Though there are no hanging
735 * node constraints in DG discretization, we define an empty
736 * AffineConstraints object that allows us to use the
737 * AffineConstraints::distribute_local_to_global() functionality.
738 *
739 * @code
740 *   AffineConstraints<double> constraints;
741 *   constraints.close();
742 *   const auto copier = [&](const CopyData &c) {
743 *   constraints.distribute_local_to_global(c.cell_matrix,
744 *   c.cell_rhs,
745 *   c.local_dof_indices,
746 *   system_matrix,
747 *   system_rhs);
748 *  
749 * @endcode
750 *
751 * Copy data from interior face assembly to the global matrix.
752 *
753 * @code
754 *   for (const CopyDataFace &cdf : c.face_data)
755 *   {
756 *   constraints.distribute_local_to_global(cdf.cell_matrix,
757 *   cdf.joint_dof_indices,
758 *   system_matrix);
759 *   }
760 *   };
761 *  
762 *  
763 * @endcode
764 *
765 * With the assembly functions defined, we can now create
766 * ScratchData and CopyData objects, and pass them together with
767 * the lambda functions above to MeshWorker::mesh_loop(). In
768 * addition, we need to specify that we want to assemble on
769 * interior faces exactly once.
770 *
771 * @code
772 *   const UpdateFlags cell_flags = update_values | update_gradients |
773 *   update_quadrature_points | update_JxW_values;
774 *   const UpdateFlags face_flags = update_values | update_gradients |
775 *   update_quadrature_points |
776 *   update_normal_vectors | update_JxW_values;
777 *  
778 *   ScratchData scratch_data(
779 *   mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
780 *   CopyData copy_data;
781 *  
782 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
783 *   dof_handler.end(),
784 *   cell_worker,
785 *   copier,
786 *   scratch_data,
787 *   copy_data,
788 *   MeshWorker::assemble_own_cells |
789 *   MeshWorker::assemble_boundary_faces |
790 *   MeshWorker::assemble_own_interior_faces_once,
791 *   boundary_worker,
792 *   face_worker);
793 *   }
794 *  
795 *  
796 *  
797 * @endcode
798 *
799 *
800 * <a name="step_74-Thesolveandoutput_resultsfunction"></a>
801 * <h3>The solve() and output_results() function</h3>
802 * The following two functions are entirely standard and without difficulty.
803 *
804 * @code
805 *   template <int dim>
806 *   void SIPGLaplace<dim>::solve()
807 *   {
808 *   SparseDirectUMFPACK A_direct;
809 *   A_direct.initialize(system_matrix);
810 *   A_direct.vmult(solution, system_rhs);
811 *   }
812 *  
813 *  
814 *  
815 *   template <int dim>
816 *   void SIPGLaplace<dim>::output_results(const unsigned int cycle) const
817 *   {
818 *   const std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) +
819 *   "-" + Utilities::int_to_string(cycle, 2) +
820 *   ".vtu";
821 *   std::ofstream output(filename);
822 *  
823 *   DataOut<dim> data_out;
824 *   data_out.attach_dof_handler(dof_handler);
825 *   data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
826 *   data_out.build_patches(mapping);
827 *   data_out.write_vtu(output);
828 *   }
829 *  
830 *  
831 * @endcode
832 *
833 *
834 * <a name="step_74-Thecompute_error_estimatefunction"></a>
835 * <h3>The compute_error_estimate() function</h3>
836 * The assembly of the error estimator here is quite similar to
837 * that of the global matrix and right-had side and can be handled
838 * by the MeshWorker::mesh_loop() framework. To understand what
839 * each of the local (lambda) functions is doing, recall first that
840 * the local cell residual is defined as
841 * @f$h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2@f$:
842 *
843 * @code
844 *   template <int dim>
845 *   void SIPGLaplace<dim>::compute_error_estimate()
846 *   {
847 *   const auto cell_worker =
848 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
849 *   ScratchData &scratch_data,
850 *   CopyData &copy_data) {
851 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
852 *  
853 *   copy_data.cell_index = cell->active_cell_index();
854 *  
855 *   const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
856 *   const unsigned int n_q_points = q_points.size();
857 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
858 *  
859 *   std::vector<Tensor<2, dim>> hessians(n_q_points);
860 *   fe_v.get_function_hessians(solution, hessians);
861 *  
862 *   std::vector<double> rhs(n_q_points);
863 *   rhs_function->value_list(q_points, rhs);
864 *  
865 *   const double hk = cell->diameter();
866 *   double residual_norm_square = 0;
867 *  
868 *   for (unsigned int point = 0; point < n_q_points; ++point)
869 *   {
870 *   const double residual =
871 *   rhs[point] + diffusion_coefficient * trace(hessians[point]);
872 *   residual_norm_square += residual * residual * JxW[point];
873 *   }
874 *   copy_data.value = hk * hk * residual_norm_square;
875 *   };
876 *  
877 * @endcode
878 *
879 * Next compute boundary terms @f$\sum_{f\in \partial K \cap \partial \Omega}
880 * \sigma \left\| [ u_h-g_D ] \right\|_f^2 @f$:
881 *
882 * @code
883 *   const auto boundary_worker =
884 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
885 *   const unsigned int &face_no,
886 *   ScratchData &scratch_data,
887 *   CopyData &copy_data) {
888 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
889 *  
890 *   const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
891 *   const unsigned n_q_points = q_points.size();
892 *  
893 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
894 *  
895 *   std::vector<double> g(n_q_points);
896 *   exact_solution->value_list(q_points, g);
897 *  
898 *   std::vector<double> sol_u(n_q_points);
899 *   fe_fv.get_function_values(solution, sol_u);
900 *  
901 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
902 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
903 *  
904 *   double difference_norm_square = 0.;
905 *   for (unsigned int point = 0; point < q_points.size(); ++point)
906 *   {
907 *   const double diff = (g[point] - sol_u[point]);
908 *   difference_norm_square += diff * diff * JxW[point];
909 *   }
910 *   copy_data.value += penalty * difference_norm_square;
911 *   };
912 *  
913 * @endcode
914 *
915 * And finally interior face terms @f$\sum_{f\in \partial K}\lbrace \sigma
916 * \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot
917 * \mathbf n ] \right\|_f^2 \rbrace@f$:
918 *
919 * @code
920 *   const auto face_worker =
921 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
922 *   const unsigned int &f,
923 *   const unsigned int &sf,
924 *   const typename DoFHandler<dim>::cell_iterator &ncell,
925 *   const unsigned int &nf,
926 *   const unsigned int &nsf,
927 *   ScratchData &scratch_data,
928 *   CopyData &copy_data) {
929 *   const FEInterfaceValues<dim> &fe_iv =
930 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
931 *  
932 *   copy_data.face_data.emplace_back();
933 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
934 *  
935 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
936 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
937 *  
938 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
939 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
940 *  
941 *   const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
942 *   const unsigned int n_q_points = q_points.size();
943 *  
944 *   std::vector<double> jump(n_q_points);
945 *   fe_iv.get_jump_in_function_values(solution, jump);
946 *  
947 *   std::vector<Tensor<1, dim>> grad_jump(n_q_points);
948 *   fe_iv.get_jump_in_function_gradients(solution, grad_jump);
949 *  
950 *   const double h = cell->face(f)->diameter();
951 *  
952 *   const double extent1 = cell->measure() / cell->face(f)->measure();
953 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
954 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
955 *  
956 *   double flux_jump_square = 0;
957 *   double u_jump_square = 0;
958 *   for (unsigned int point = 0; point < n_q_points; ++point)
959 *   {
960 *   u_jump_square += jump[point] * jump[point] * JxW[point];
961 *   const double flux_jump = grad_jump[point] * normals[point];
962 *   flux_jump_square +=
963 *   diffusion_coefficient * flux_jump * flux_jump * JxW[point];
964 *   }
965 *   copy_data_face.values[0] =
966 *   0.5 * h * (flux_jump_square + penalty * u_jump_square);
967 *   copy_data_face.values[1] = copy_data_face.values[0];
968 *   };
969 *  
970 * @endcode
971 *
972 * Having computed local contributions for each cell, we still
973 * need a way to copy these into the global vector that will hold
974 * the error estimators for all cells:
975 *
976 * @code
977 *   const auto copier = [&](const CopyData &copy_data) {
978 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
979 *   estimated_error_square_per_cell[copy_data.cell_index] +=
980 *   copy_data.value;
981 *   for (const CopyDataFace &cdf : copy_data.face_data)
982 *   for (unsigned int j = 0; j < 2; ++j)
983 *   estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
984 *   };
985 *  
986 * @endcode
987 *
988 * After all of this set-up, let's do the actual work: We resize
989 * the vector into which the results will be written, and then
990 * drive the whole process using the MeshWorker::mesh_loop()
991 * function.
992 *
993 * @code
994 *   estimated_error_square_per_cell.reinit(triangulation.n_active_cells());
995 *  
996 *   const UpdateFlags cell_flags =
998 *   const UpdateFlags face_flags = update_values | update_gradients |
1001 *  
1002 *   ScratchData scratch_data(
1003 *   mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
1004 *  
1005 *   CopyData copy_data;
1006 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1007 *   dof_handler.end(),
1008 *   cell_worker,
1009 *   copier,
1010 *   scratch_data,
1011 *   copy_data,
1015 *   boundary_worker,
1016 *   face_worker);
1017 *   }
1018 *  
1019 * @endcode
1020 *
1021 *
1022 * <a name="step_74-Thecompute_energy_norm_errorfunction"></a>
1023 * <h3>The compute_energy_norm_error() function</h3>
1024 * Next, we evaluate the accuracy in terms of the energy norm.
1025 * This function is similar to the assembling of the error estimator above.
1026 * Here we compute the square of the energy norm defined by
1027 * @f[
1028 * \|u \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla u \|_K^2 +
1029 * \sum_{f \in F_i} \sigma \| [ u ] \|_f^2 +
1030 * \sum_{f \in F_b} \sigma \|u\|_f^2.
1031 * @f]
1032 * Therefore the corresponding error is
1033 * @f[
1034 * \|u -u_h \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2
1035 * + \sum_{f \in F_i} \sigma \|[ u_h ] \|_f^2 + \sum_{f \in F_b}\sigma
1036 * \|u_h-g_D\|_f^2.
1037 * @f]
1038 *
1039 * @code
1040 *   template <int dim>
1041 *   double SIPGLaplace<dim>::compute_energy_norm_error()
1042 *   {
1043 *   energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
1044 *  
1045 * @endcode
1046 *
1047 * Assemble @f$\sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2 @f$.
1048 *
1049 * @code
1050 *   const auto cell_worker =
1051 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
1052 *   ScratchData &scratch_data,
1053 *   CopyData &copy_data) {
1054 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
1055 *  
1056 *   copy_data.cell_index = cell->active_cell_index();
1057 *  
1058 *   const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
1059 *   const unsigned int n_q_points = q_points.size();
1060 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
1061 *  
1062 *   std::vector<Tensor<1, dim>> grad_u(n_q_points);
1063 *   fe_v.get_function_gradients(solution, grad_u);
1064 *  
1065 *   std::vector<Tensor<1, dim>> grad_exact(n_q_points);
1066 *   exact_solution->gradient_list(q_points, grad_exact);
1067 *  
1068 *   double norm_square = 0;
1069 *   for (unsigned int point = 0; point < n_q_points; ++point)
1070 *   {
1071 *   norm_square +=
1072 *   (grad_u[point] - grad_exact[point]).norm_square() * JxW[point];
1073 *   }
1074 *   copy_data.value = diffusion_coefficient * norm_square;
1075 *   };
1076 *  
1077 * @endcode
1078 *
1079 * Assemble @f$\sum_{f \in F_b}\sigma \|u_h-g_D\|_f^2@f$.
1080 *
1081 * @code
1082 *   const auto boundary_worker =
1083 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
1084 *   const unsigned int &face_no,
1085 *   ScratchData &scratch_data,
1086 *   CopyData &copy_data) {
1087 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
1088 *  
1089 *   const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
1090 *   const unsigned n_q_points = q_points.size();
1091 *  
1092 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
1093 *  
1094 *   std::vector<double> g(n_q_points);
1095 *   exact_solution->value_list(q_points, g);
1096 *  
1097 *   std::vector<double> sol_u(n_q_points);
1098 *   fe_fv.get_function_values(solution, sol_u);
1099 *  
1100 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
1101 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
1102 *  
1103 *   double difference_norm_square = 0.;
1104 *   for (unsigned int point = 0; point < q_points.size(); ++point)
1105 *   {
1106 *   const double diff = (g[point] - sol_u[point]);
1107 *   difference_norm_square += diff * diff * JxW[point];
1108 *   }
1109 *   copy_data.value += penalty * difference_norm_square;
1110 *   };
1111 *  
1112 * @endcode
1113 *
1114 * Assemble @f$\sum_{f \in F_i} \sigma \| [ u_h ] \|_f^2@f$.
1115 *
1116 * @code
1117 *   const auto face_worker =
1118 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
1119 *   const unsigned int &f,
1120 *   const unsigned int &sf,
1121 *   const typename DoFHandler<dim>::cell_iterator &ncell,
1122 *   const unsigned int &nf,
1123 *   const unsigned int &nsf,
1124 *   ScratchData &scratch_data,
1125 *   CopyData &copy_data) {
1126 *   const FEInterfaceValues<dim> &fe_iv =
1127 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
1128 *  
1129 *   copy_data.face_data.emplace_back();
1130 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
1131 *  
1132 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
1133 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
1134 *  
1135 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
1136 *  
1137 *   const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
1138 *   const unsigned int n_q_points = q_points.size();
1139 *  
1140 *   std::vector<double> jump(n_q_points);
1141 *   fe_iv.get_jump_in_function_values(solution, jump);
1142 *  
1143 *   const double extent1 = cell->measure() / cell->face(f)->measure();
1144 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
1145 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
1146 *  
1147 *   double u_jump_square = 0;
1148 *   for (unsigned int point = 0; point < n_q_points; ++point)
1149 *   {
1150 *   u_jump_square += jump[point] * jump[point] * JxW[point];
1151 *   }
1152 *   copy_data_face.values[0] = 0.5 * penalty * u_jump_square;
1153 *   copy_data_face.values[1] = copy_data_face.values[0];
1154 *   };
1155 *  
1156 *   const auto copier = [&](const CopyData &copy_data) {
1157 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
1158 *   energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
1159 *   for (const CopyDataFace &cdf : copy_data.face_data)
1160 *   for (unsigned int j = 0; j < 2; ++j)
1161 *   energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1162 *   };
1163 *  
1164 *   const UpdateFlags cell_flags =
1166 *   UpdateFlags face_flags =
1168 *  
1169 *   const ScratchData scratch_data(mapping,
1170 *   fe,
1171 *   quadrature_overintegration,
1172 *   cell_flags,
1173 *   face_quadrature_overintegration,
1174 *   face_flags);
1175 *  
1176 *   CopyData copy_data;
1177 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1178 *   dof_handler.end(),
1179 *   cell_worker,
1180 *   copier,
1181 *   scratch_data,
1182 *   copy_data,
1186 *   boundary_worker,
1187 *   face_worker);
1188 *   const double energy_error =
1189 *   std::sqrt(energy_norm_square_per_cell.l1_norm());
1190 *   return energy_error;
1191 *   }
1192 *  
1193 *  
1194 *  
1195 * @endcode
1196 *
1197 *
1198 * <a name="step_74-Therefine_gridfunction"></a>
1199 * <h3>The refine_grid() function</h3>
1200 *
1201 * @code
1202 *   template <int dim>
1203 *   void SIPGLaplace<dim>::refine_grid()
1204 *   {
1205 *   const double refinement_fraction = 0.1;
1206 *  
1208 *   triangulation, estimated_error_square_per_cell, refinement_fraction, 0.);
1209 *  
1211 *   }
1212 *  
1213 *  
1214 *  
1215 * @endcode
1216 *
1217 *
1218 * <a name="step_74-Thecompute_errorsfunction"></a>
1219 * <h3>The compute_errors() function</h3>
1220 * We compute three errors in the @f$L_2@f$ norm, @f$H_1@f$ seminorm, and
1221 * the energy norm, respectively. These are then printed to screen,
1222 * but also stored in a table that records how these errors decay
1223 * with mesh refinement and which can be output in one step at the
1224 * end of the program.
1225 *
1226 * @code
1227 *   template <int dim>
1228 *   void SIPGLaplace<dim>::compute_errors()
1229 *   {
1230 *   double L2_error, H1_error, energy_error;
1231 *  
1232 *   {
1233 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
1235 *   dof_handler,
1236 *   solution,
1237 *   *(exact_solution.get()),
1238 *   difference_per_cell,
1239 *   quadrature_overintegration,
1241 *  
1243 *   difference_per_cell,
1245 *   convergence_table.add_value("L2", L2_error);
1246 *   }
1247 *  
1248 *   {
1249 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
1251 *   dof_handler,
1252 *   solution,
1253 *   *(exact_solution.get()),
1254 *   difference_per_cell,
1255 *   quadrature_overintegration,
1257 *  
1259 *   difference_per_cell,
1261 *   convergence_table.add_value("H1", H1_error);
1262 *   }
1263 *  
1264 *   {
1265 *   energy_error = compute_energy_norm_error();
1266 *   convergence_table.add_value("Energy", energy_error);
1267 *   }
1268 *  
1269 *   std::cout << " Error in the L2 norm : " << L2_error << std::endl
1270 *   << " Error in the H1 seminorm : " << H1_error << std::endl
1271 *   << " Error in the energy norm : " << energy_error
1272 *   << std::endl;
1273 *   }
1274 *  
1275 *  
1276 *  
1277 * @endcode
1278 *
1279 *
1280 * <a name="step_74-Therunfunction"></a>
1281 * <h3>The run() function</h3>
1282 *
1283 * @code
1284 *   template <int dim>
1285 *   void SIPGLaplace<dim>::run()
1286 *   {
1287 *   const unsigned int max_cycle =
1288 *   (test_case == TestCase::convergence_rate ? 6 : 20);
1289 *   for (unsigned int cycle = 0; cycle < max_cycle; ++cycle)
1290 *   {
1291 *   std::cout << "Cycle " << cycle << std::endl;
1292 *  
1293 *   switch (test_case)
1294 *   {
1295 *   case TestCase::convergence_rate:
1296 *   {
1297 *   if (cycle == 0)
1298 *   {
1300 *  
1302 *   }
1303 *   else
1304 *   {
1306 *   }
1307 *   break;
1308 *   }
1309 *  
1310 *   case TestCase::l_singularity:
1311 *   {
1312 *   if (cycle == 0)
1313 *   {
1316 *   }
1317 *   else
1318 *   {
1319 *   refine_grid();
1320 *   }
1321 *   break;
1322 *   }
1323 *  
1324 *   default:
1325 *   {
1327 *   }
1328 *   }
1329 *  
1330 *   std::cout << " Number of active cells : "
1331 *   << triangulation.n_active_cells() << std::endl;
1332 *   setup_system();
1333 *  
1334 *   std::cout << " Number of degrees of freedom : " << dof_handler.n_dofs()
1335 *   << std::endl;
1336 *  
1337 *   assemble_system();
1338 *   solve();
1339 *   output_results(cycle);
1340 *   {
1341 *   convergence_table.add_value("cycle", cycle);
1342 *   convergence_table.add_value("cells", triangulation.n_active_cells());
1343 *   convergence_table.add_value("dofs", dof_handler.n_dofs());
1344 *   }
1345 *   compute_errors();
1346 *  
1347 *   if (test_case == TestCase::l_singularity)
1348 *   {
1349 *   compute_error_estimate();
1350 *   std::cout << " Estimated error : "
1351 *   << std::sqrt(estimated_error_square_per_cell.l1_norm())
1352 *   << std::endl;
1353 *  
1354 *   convergence_table.add_value(
1355 *   "Estimator",
1356 *   std::sqrt(estimated_error_square_per_cell.l1_norm()));
1357 *   }
1358 *   std::cout << std::endl;
1359 *   }
1360 *  
1361 * @endcode
1362 *
1363 * Having run all of our computations, let us tell the convergence
1364 * table how to format its data and output it to screen:
1365 *
1366 * @code
1367 *   convergence_table.set_precision("L2", 3);
1368 *   convergence_table.set_precision("H1", 3);
1369 *   convergence_table.set_precision("Energy", 3);
1370 *  
1371 *   convergence_table.set_scientific("L2", true);
1372 *   convergence_table.set_scientific("H1", true);
1373 *   convergence_table.set_scientific("Energy", true);
1374 *  
1375 *   if (test_case == TestCase::convergence_rate)
1376 *   {
1377 *   convergence_table.evaluate_convergence_rates(
1379 *   convergence_table.evaluate_convergence_rates(
1381 *   }
1382 *   if (test_case == TestCase::l_singularity)
1383 *   {
1384 *   convergence_table.set_precision("Estimator", 3);
1385 *   convergence_table.set_scientific("Estimator", true);
1386 *   }
1387 *  
1388 *   std::cout << "degree = " << degree << std::endl;
1389 *   convergence_table.write_text(
1391 *   }
1392 *   } // namespace Step74
1393 *  
1394 *  
1395 *  
1396 * @endcode
1397 *
1398 *
1399 * <a name="step_74-Themainfunction"></a>
1400 * <h3>The main() function</h3>
1401 * The following <code>main</code> function is similar to previous examples as
1402 * well, and need not be commented on.
1403 *
1404 * @code
1405 *   int main()
1406 *   {
1407 *   try
1408 *   {
1409 *   using namespace dealii;
1410 *   using namespace Step74;
1411 *  
1412 *   const TestCase test_case = TestCase::l_singularity;
1413 *  
1414 *   SIPGLaplace<2> problem(test_case);
1415 *   problem.run();
1416 *   }
1417 *   catch (std::exception &exc)
1418 *   {
1419 *   std::cerr << std::endl
1420 *   << std::endl
1421 *   << "----------------------------------------------------"
1422 *   << std::endl;
1423 *   std::cerr << "Exception on processing: " << std::endl
1424 *   << exc.what() << std::endl
1425 *   << "Aborting!" << std::endl
1426 *   << "----------------------------------------------------"
1427 *   << std::endl;
1428 *   return 1;
1429 *   }
1430 *   catch (...)
1431 *   {
1432 *   std::cerr << std::endl
1433 *   << std::endl
1434 *   << "----------------------------------------------------"
1435 *   << std::endl;
1436 *   std::cerr << "Unknown exception!" << std::endl
1437 *   << "Aborting!" << std::endl
1438 *   << "----------------------------------------------------"
1439 *   << std::endl;
1440 *   return 1;
1441 *   };
1442 *  
1443 *   return 0;
1444 *   }
1445 * @endcode
1446<a name="step_74-Results"></a><h1>Results</h1>
1447
1448
1449The output of this program consist of the console output and
1450solutions in vtu format.
1451
1452In the first test case, when you run the program, the screen output should look like the following:
1453@code
1454Cycle 0
1455 Number of active cells : 16
1456 Number of degrees of freedom : 256
1457 Error in the L2 norm : 0.00193285
1458 Error in the H1 seminorm : 0.106087
1459 Error in the energy norm : 0.150625
1460
1461Cycle 1
1462 Number of active cells : 64
1463 Number of degrees of freedom : 1024
1464 Error in the L2 norm : 9.60497e-05
1465 Error in the H1 seminorm : 0.0089954
1466 Error in the energy norm : 0.0113265
1467
1468Cycle 2
1469.
1470.
1471.
1472@endcode
1473
1474When using the smooth case with polynomial degree 3, the convergence
1475table will look like this:
1476<table align="center" class="doxtable">
1477 <tr>
1478 <th>cycle</th>
1479 <th>n_cells</th>
1480 <th>n_dofs</th>
1481 <th>L2 </th>
1482 <th>rate</th>
1483 <th>H1</th>
1484 <th>rate</th>
1485 <th>Energy</th>
1486 </tr>
1487 <tr>
1488 <td align="center">0</td>
1489 <td align="right">16</td>
1490 <td align="right">256</td>
1491 <td align="center">1.933e-03</td>
1492 <td>&nbsp;</td>
1493 <td align="center">1.061e-01</td>
1494 <td>&nbsp;</td>
1495 <td align="center">1.506e-01</td>
1496 </tr>
1497 <tr>
1498 <td align="center">1</td>
1499 <td align="right">64</td>
1500 <td align="right">1024</td>
1501 <td align="center">9.605e-05</td>
1502 <td align="center">4.33</td>
1503 <td align="center">8.995e-03</td>
1504 <td align="center">3.56</td>
1505 <td align="center">1.133e-02</td>
1506 </tr>
1507 <tr>
1508 <td align="center">2</td>
1509 <td align="right">256</td>
1510 <td align="right">4096</td>
1511 <td align="center">5.606e-06</td>
1512 <td align="center">4.10</td>
1513 <td align="center">9.018e-04</td>
1514 <td align="center">3.32</td>
1515 <td align="center">9.736e-04</td>
1516 </tr>
1517 <tr>
1518 <td align="center">3</td>
1519 <td align="right">1024</td>
1520 <td align="right">16384</td>
1521 <td align="center">3.484e-07</td>
1522 <td align="center">4.01</td>
1523 <td align="center">1.071e-04</td>
1524 <td align="center">3.07</td>
1525 <td align="center">1.088e-04</td>
1526 </tr>
1527 <tr>
1528 <td align="center">4</td>
1529 <td align="right">4096</td>
1530 <td align="right">65536</td>
1531 <td align="center">2.179e-08</td>
1532 <td align="center">4.00</td>
1533 <td align="center">1.327e-05</td>
1534 <td align="center">3.01</td>
1535 <td align="center">1.331e-05</td>
1536 </tr>
1537 <tr>
1538 <td align="center">5</td>
1539 <td align="right">16384</td>
1540 <td align="right">262144</td>
1541 <td align="center">1.363e-09</td>
1542 <td align="center">4.00</td>
1543 <td align="center">1.656e-06</td>
1544 <td align="center">3.00</td>
1545 <td align="center">1.657e-06</td>
1546 </tr>
1547</table>
1548
1549Theoretically, for polynomial degree @f$p@f$, the order of convergence in @f$L_2@f$
1550norm and @f$H^1@f$ seminorm should be @f$p+1@f$ and @f$p@f$, respectively. Our numerical
1551results are in good agreement with theory.
1552
1553In the second test case, when you run the program, the screen output should look like the following:
1554@code
1555Cycle 0
1556 Number of active cells : 192
1557 Number of degrees of freedom : 3072
1558 Error in the L2 norm : 0.000323585
1559 Error in the H1 seminorm : 0.0296202
1560 Error in the energy norm : 0.0420478
1561 Estimated error : 0.136067
1562
1563Cycle 1
1564 Number of active cells : 249
1565 Number of degrees of freedom : 3984
1566 Error in the L2 norm : 0.000114739
1567 Error in the H1 seminorm : 0.0186571
1568 Error in the energy norm : 0.0264879
1569 Estimated error : 0.0857186
1570
1571Cycle 2
1572.
1573.
1574.
1575@endcode
1576
1577The following figure provides a log-log plot of the errors versus
1578the number of degrees of freedom for this test case on the L-shaped
1579domain. In order to interpret it, let @f$n@f$ be the number of degrees of
1580freedom, then on uniformly refined meshes, @f$h@f$ is of order
1581@f$1/\sqrt{n}@f$ in 2D. Combining the theoretical results in the previous case,
1582we see that if the solution is sufficiently smooth,
1583we can expect the error in the @f$L_2@f$ norm to be of order @f$O(n^{-\frac{p+1}{2}})@f$
1584and in @f$H^1@f$ seminorm to be @f$O(n^{-\frac{p}{2}})@f$. It is not a priori
1585clear that one would get the same kind of behavior as a function of
1586@f$n@f$ on adaptively refined meshes like the ones we use for this second
1587test case, but one can certainly hope. Indeed, from the figure, we see
1588that the SIPG with adaptive mesh refinement produces asymptotically
1589the kinds of hoped-for results:
1590
1591<img width="600px" src="https://www.dealii.org/images/steps/developer/step-74.log-log-plot.png" alt="">
1592
1593In addition, we observe that the error estimator decreases
1594at almost the same rate as the errors in the energy norm and @f$H^1@f$ seminorm,
1595and one order lower than the @f$L_2@f$ error. This suggests
1596its ability to predict regions with large errors.
1597
1598While this tutorial is focused on the implementation, the @ref step_59 "step-59" tutorial program achieves an efficient
1599large-scale solver in terms of computing time with matrix-free solution techniques.
1600Note that the @ref step_59 "step-59" tutorial does not work with meshes containing hanging nodes at this moment,
1601because the multigrid interface matrices are not as easily determined,
1602but that is merely the lack of some interfaces in deal.II, nothing fundamental.
1603 *
1604 *
1605<a name="step_74-PlainProg"></a>
1606<h1> The plain program</h1>
1607@include "step-74.cc"
1608*/
void reinit(const CellIteratorType &cell, const unsigned int face_no, const unsigned int sub_face_no, const CellNeighborIteratorType &cell_neighbor, const unsigned int face_no_neighbor, const unsigned int sub_face_no_neighbor, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int, const unsigned int fe_index_neighbor=numbers::invalid_unsigned_int)
const std::vector< Point< spacedim > > & get_quadrature_points() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
Definition point.h:111
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int cell_index
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:281
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition l2.h:159
@ assemble_boundary_faces
@ assemble_own_interior_faces_once
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
void free(T *&pointer)
Definition cuda.h:96
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14883
int(& functions)(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
Definition numbers.h:259
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation