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