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-82.h
Go to the documentation of this file.
1
524) const
525 *   {
526 *   double return_value = 0.0;
527 *  
528 *   if (dim == 2)
529 *   {
530 *   return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
531 *   +24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0])) +
532 *   2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
533 *   (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]);
534 *   }
535 *   else if (dim == 3)
536 *   {
537 *   return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) *
538 *   p[2] * (1.0 - p[2])) +
539 *   24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
540 *   p[2] * (1.0 - p[2])) +
541 *   24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
542 *   p[1] * (1.0 - p[1])) +
543 *   2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
544 *   (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
545 *   Utilities::fixed_power<2>(p[2] * (1.0 - p[2])) +
546 *   2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
547 *   (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
548 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
549 *   2.0 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
550 *   (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
551 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
552 *   }
553 *   else
555 *  
556 *   return return_value;
557 *   }
558 *  
559 *  
560 *  
561 * @endcode
562 *
563 * This class implement the manufactured (exact) solution @f$u@f$. To compute the
564 * errors, we need the value of @f$u@f$ as well as its gradient and its Hessian.
565 *
566 * @code
567 *   template <int dim>
568 *   class ExactSolution : public Function<dim>
569 *   {
570 *   public:
571 *   ExactSolution()
572 *   : Function<dim>()
573 *   {}
574 *  
575 *   virtual double value(const Point<dim> &p,
576 *   const unsigned int component = 0) const override;
577 *  
578 *   virtual Tensor<1, dim>
579 *   gradient(const Point<dim> &p,
580 *   const unsigned int component = 0) const override;
581 *  
582 *   virtual SymmetricTensor<2, dim>
583 *   hessian(const Point<dim> &p,
584 *   const unsigned int component = 0) const override;
585 *   };
586 *  
587 *  
588 *  
589 *   template <int dim>
590 *   double ExactSolution<dim>::value(const Point<dim> &p,
591 *   const unsigned int /*component*/) const
592 *   {
593 *   double return_value = 0.0;
594 *  
595 *   if (dim == 2)
596 *   {
597 *   return_value =
598 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
599 *   }
600 *   else if (dim == 3)
601 *   {
602 *   return_value = Utilities::fixed_power<2>(
603 *   p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
604 *   }
605 *   else
607 *  
608 *   return return_value;
609 *   }
610 *  
611 *  
612 *  
613 *   template <int dim>
614 *   Tensor<1, dim>
615 *   ExactSolution<dim>::gradient(const Point<dim> &p,
616 *   const unsigned int /*component*/) const
617 *   {
618 *   Tensor<1, dim> return_gradient;
619 *  
620 *   if (dim == 2)
621 *   {
622 *   return_gradient[0] =
623 *   (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
624 *   4.0 * Utilities::fixed_power<3>(p[0])) *
625 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
626 *   return_gradient[1] =
627 *   (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
628 *   4.0 * Utilities::fixed_power<3>(p[1])) *
629 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
630 *   }
631 *   else if (dim == 3)
632 *   {
633 *   return_gradient[0] =
634 *   (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
635 *   4.0 * Utilities::fixed_power<3>(p[0])) *
636 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
637 *   return_gradient[1] =
638 *   (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
639 *   4.0 * Utilities::fixed_power<3>(p[1])) *
640 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
641 *   return_gradient[2] =
642 *   (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
643 *   4.0 * Utilities::fixed_power<3>(p[2])) *
644 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
645 *   }
646 *   else
648 *  
649 *   return return_gradient;
650 *   }
651 *  
652 *  
653 *  
654 *   template <int dim>
656 *   ExactSolution<dim>::hessian(const Point<dim> &p,
657 *   const unsigned int /*component*/) const
658 *   {
659 *   SymmetricTensor<2, dim> return_hessian;
660 *  
661 *   if (dim == 2)
662 *   {
663 *   return_hessian[0][0] = (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
664 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
665 *   return_hessian[0][1] =
666 *   (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
667 *   4.0 * Utilities::fixed_power<3>(p[0])) *
668 *   (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
669 *   4.0 * Utilities::fixed_power<3>(p[1]));
670 *   return_hessian[1][1] = (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
671 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
672 *   }
673 *   else if (dim == 3)
674 *   {
675 *   return_hessian[0][0] =
676 *   (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
677 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
678 *   return_hessian[0][1] =
679 *   (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
680 *   4.0 * Utilities::fixed_power<3>(p[0])) *
681 *   (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
682 *   4.0 * Utilities::fixed_power<3>(p[1])) *
683 *   Utilities::fixed_power<2>(p[2] * (1.0 - p[2]));
684 *   return_hessian[0][2] =
685 *   (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
686 *   4.0 * Utilities::fixed_power<3>(p[0])) *
687 *   (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
688 *   4.0 * Utilities::fixed_power<3>(p[2])) *
689 *   Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
690 *   return_hessian[1][1] =
691 *   (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
692 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
693 *   return_hessian[1][2] =
694 *   (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
695 *   4.0 * Utilities::fixed_power<3>(p[1])) *
696 *   (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
697 *   4.0 * Utilities::fixed_power<3>(p[2])) *
698 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
699 *   return_hessian[2][2] =
700 *   (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
701 *   Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
702 *   }
703 *   else
705 *  
706 *   return return_hessian;
707 *   }
708 *  
709 *  
710 *  
711 * @endcode
712 *
713 *
714 * <a name="step_82-ImplementationofthecodeBiLaplacianLDGLiftcodeclass"></a>
715 * <h3>Implementation of the <code>BiLaplacianLDGLift</code> class</h3>
716 *
717
718 *
719 *
720 * <a name="step_82-BiLaplacianLDGLiftBiLaplacianLDGLift"></a>
721 * <h4>BiLaplacianLDGLift::BiLaplacianLDGLift</h4>
722 *
723
724 *
725 * In the constructor, we set the polynomial degree of the two finite element
726 * spaces, we associate the corresponding DoF handlers to the triangulation,
727 * and we set the two penalty coefficients.
728 *
729 * @code
730 *   template <int dim>
731 *   BiLaplacianLDGLift<dim>::BiLaplacianLDGLift(const unsigned int n_refinements,
732 *   const unsigned int fe_degree,
733 *   const double penalty_jump_grad,
734 *   const double penalty_jump_val)
735 *   : n_refinements(n_refinements)
736 *   , fe(fe_degree)
737 *   , dof_handler(triangulation)
738 *   , fe_lift(FE_DGQ<dim>(fe_degree), dim * dim)
739 *   , penalty_jump_grad(penalty_jump_grad)
740 *   , penalty_jump_val(penalty_jump_val)
741 *   {}
742 *  
743 *  
744 *  
745 * @endcode
746 *
747 *
748 * <a name="step_82-BiLaplacianLDGLiftmake_grid"></a>
749 * <h4>BiLaplacianLDGLift::make_grid</h4>
750 *
751
752 *
753 * To build a mesh for @f$\Omega=(0,1)^d@f$, we simply call the function
754 * <code>GridGenerator::hyper_cube</code> and then refine it using
755 * <code>refine_global</code>. The number of refinements is hard-coded
756 * here.
757 *
758 * @code
759 *   template <int dim>
760 *   void BiLaplacianLDGLift<dim>::make_grid()
761 *   {
762 *   std::cout << "Building the mesh............." << std::endl;
763 *  
765 *  
766 *   triangulation.refine_global(n_refinements);
767 *  
768 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
769 *   << std::endl;
770 *   }
771 *  
772 *  
773 *  
774 * @endcode
775 *
776 *
777 * <a name="step_82-BiLaplacianLDGLiftsetup_system"></a>
778 * <h4>BiLaplacianLDGLift::setup_system</h4>
779 *
780
781 *
782 * In the following function, we set up the degrees of freedom, the sparsity
783 * pattern, the size of the matrix @f$A@f$, and the size of the solution and
784 * right-hand side vectors
785 * @f$\boldsymbol{U}@f$ and @f$\boldsymbol{F}@f$. For the sparsity pattern, we cannot
786 * directly use the function <code>DoFTools::make_flux_sparsity_pattern</code>
787 * (as we would do for instance for the SIPG method) because we need to take
788 * into account the interactions of a neighboring cell with another
789 * neighboring cell as described in the introduction. The extended sparsity
790 * pattern is built by iterating over all the active cells. For the current
791 * cell, we collect all its degrees of freedom as well as the degrees of
792 * freedom of all its neighboring cells, and then couple everything with
793 * everything.
794 *
795 * @code
796 *   template <int dim>
797 *   void BiLaplacianLDGLift<dim>::setup_system()
798 *   {
799 *   dof_handler.distribute_dofs(fe);
800 *  
801 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
802 *   << std::endl;
803 *  
804 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
805 *  
806 *   const auto dofs_per_cell = fe.dofs_per_cell;
807 *  
808 *   for (const auto &cell : dof_handler.active_cell_iterators())
809 *   {
810 *   std::vector<types::global_dof_index> dofs(dofs_per_cell);
811 *   cell->get_dof_indices(dofs);
812 *  
813 *   for (unsigned int f = 0; f < cell->n_faces(); ++f)
814 *   if (!cell->face(f)->at_boundary())
815 *   {
816 *   const auto neighbor_cell = cell->neighbor(f);
817 *  
818 *   std::vector<types::global_dof_index> tmp(dofs_per_cell);
819 *   neighbor_cell->get_dof_indices(tmp);
820 *  
821 *   dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp));
822 *   }
823 *  
824 *   for (const auto i : dofs)
825 *   for (const auto j : dofs)
826 *   {
827 *   dsp.add(i, j);
828 *   dsp.add(j, i);
829 *   }
830 *   }
831 *  
832 *   sparsity_pattern.copy_from(dsp);
833 *  
834 *  
835 *   matrix.reinit(sparsity_pattern);
836 *   rhs.reinit(dof_handler.n_dofs());
837 *  
838 *   solution.reinit(dof_handler.n_dofs());
839 *  
840 * @endcode
841 *
842 * At the end of the function, we output this sparsity pattern as
843 * a scalable vector graphic. You can visualize it by loading this
844 * file in most web browsers:
845 *
846 * @code
847 *   std::ofstream out("sparsity-pattern.svg");
848 *   sparsity_pattern.print_svg(out);
849 *   }
850 *  
851 *  
852 *  
853 * @endcode
854 *
855 *
856 * <a name="step_82-BiLaplacianLDGLiftassemble_system"></a>
857 * <h4>BiLaplacianLDGLift::assemble_system</h4>
858 *
859
860 *
861 * This function simply calls the two functions responsible
862 * for the assembly of the matrix and the right-hand side.
863 *
864 * @code
865 *   template <int dim>
866 *   void BiLaplacianLDGLift<dim>::assemble_system()
867 *   {
868 *   std::cout << "Assembling the system............." << std::endl;
869 *  
870 *   assemble_matrix();
871 *   assemble_rhs();
872 *  
873 *   std::cout << "Done. " << std::endl;
874 *   }
875 *  
876 *  
877 *  
878 * @endcode
879 *
880 *
881 * <a name="step_82-BiLaplacianLDGLiftassemble_matrix"></a>
882 * <h4>BiLaplacianLDGLift::assemble_matrix</h4>
883 *
884
885 *
886 * This function assembles the matrix @f$A@f$ whose entries are defined
887 * by @f$A_{ij}=A_h(\varphi_j,\varphi_i)@f$ which involves the product of
888 * discrete Hessians and the penalty terms.
889 *
890 * @code
891 *   template <int dim>
892 *   void BiLaplacianLDGLift<dim>::assemble_matrix()
893 *   {
894 *   matrix = 0;
895 *  
896 *   const QGauss<dim> quad(fe.degree + 1);
897 *   const QGauss<dim - 1> quad_face(fe.degree + 1);
898 *  
899 *   const unsigned int n_q_points = quad.size();
900 *   const unsigned int n_q_points_face = quad_face.size();
901 *  
902 *   FEValues<dim> fe_values(fe, quad, update_hessians | update_JxW_values);
903 *  
904 *   FEFaceValues<dim> fe_face(
906 *  
907 *   FEFaceValues<dim> fe_face_neighbor(
909 *  
910 *   const unsigned int n_dofs = fe_values.dofs_per_cell;
911 *  
912 *   std::vector<types::global_dof_index> local_dof_indices(n_dofs);
913 *   std::vector<types::global_dof_index> local_dof_indices_neighbor(n_dofs);
914 *   std::vector<types::global_dof_index> local_dof_indices_neighbor_2(n_dofs);
915 *  
916 * @endcode
917 *
918 * As indicated in the introduction, the following matrices are used for
919 * the contributions of the products of the discrete Hessians.
920 *
921 * @code
922 *   FullMatrix<double> stiffness_matrix_cc(n_dofs,
923 *   n_dofs); // interactions cell / cell
924 *   FullMatrix<double> stiffness_matrix_cn(
925 *   n_dofs, n_dofs); // interactions cell / neighbor
926 *   FullMatrix<double> stiffness_matrix_nc(
927 *   n_dofs, n_dofs); // interactions neighbor / cell
928 *   FullMatrix<double> stiffness_matrix_nn(
929 *   n_dofs, n_dofs); // interactions neighbor / neighbor
930 *   FullMatrix<double> stiffness_matrix_n1n2(
931 *   n_dofs, n_dofs); // interactions neighbor1 / neighbor2
932 *   FullMatrix<double> stiffness_matrix_n2n1(
933 *   n_dofs, n_dofs); // interactions neighbor2 / neighbor1
934 *  
935 * @endcode
936 *
937 * The following matrices are used for the contributions of the two
938 * penalty terms:
939 *
940 * @code
941 *   FullMatrix<double> ip_matrix_cc(n_dofs, n_dofs); // interactions cell / cell
942 *   FullMatrix<double> ip_matrix_cn(n_dofs,
943 *   n_dofs); // interactions cell / neighbor
944 *   FullMatrix<double> ip_matrix_nc(n_dofs,
945 *   n_dofs); // interactions neighbor / cell
946 *   FullMatrix<double> ip_matrix_nn(n_dofs,
947 *   n_dofs); // interactions neighbor / neighbor
948 *  
949 *   std::vector<std::vector<Tensor<2, dim>>> discrete_hessians(
950 *   n_dofs, std::vector<Tensor<2, dim>>(n_q_points));
951 *   std::vector<std::vector<std::vector<Tensor<2, dim>>>>
952 *   discrete_hessians_neigh(GeometryInfo<dim>::faces_per_cell,
953 *   discrete_hessians);
954 *  
955 *   for (const auto &cell : dof_handler.active_cell_iterators())
956 *   {
957 *   fe_values.reinit(cell);
958 *   cell->get_dof_indices(local_dof_indices);
959 *  
960 * @endcode
961 *
962 * We now compute all the discrete Hessians that are not vanishing
963 * on the current cell, i.e., the discrete Hessian of all the basis
964 * functions with support on the current cell or on one of its
965 * neighbors.
966 *
967 * @code
968 *   compute_discrete_hessians(cell,
969 *   discrete_hessians,
970 *   discrete_hessians_neigh);
971 *  
972 * @endcode
973 *
974 * First, we compute and add the interactions of the degrees of freedom
975 * of the current cell.
976 *
977 * @code
978 *   stiffness_matrix_cc = 0;
979 *   for (unsigned int q = 0; q < n_q_points; ++q)
980 *   {
981 *   const double dx = fe_values.JxW(q);
982 *  
983 *   for (unsigned int i = 0; i < n_dofs; ++i)
984 *   for (unsigned int j = 0; j < n_dofs; ++j)
985 *   {
986 *   const Tensor<2, dim> &H_i = discrete_hessians[i][q];
987 *   const Tensor<2, dim> &H_j = discrete_hessians[j][q];
988 *  
989 *   stiffness_matrix_cc(i, j) += scalar_product(H_j, H_i) * dx;
990 *   }
991 *   }
992 *  
993 *   for (unsigned int i = 0; i < n_dofs; ++i)
994 *   for (unsigned int j = 0; j < n_dofs; ++j)
995 *   {
996 *   matrix(local_dof_indices[i], local_dof_indices[j]) +=
997 *   stiffness_matrix_cc(i, j);
998 *   }
999 *  
1000 * @endcode
1001 *
1002 * Next, we compute and add the interactions of the degrees of freedom
1003 * of the current cell with those of its neighbors. Note that the
1004 * interactions of the degrees of freedom of a neighbor with those of
1005 * the same neighbor are included here.
1006 *
1007 * @code
1008 *   for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1009 *   {
1010 *   const typename DoFHandler<dim>::face_iterator face =
1011 *   cell->face(face_no);
1012 *  
1013 *   const bool at_boundary = face->at_boundary();
1014 *   if (!at_boundary)
1015 *   {
1016 * @endcode
1017 *
1018 * There is nothing to be done if boundary face (the liftings of
1019 * the Dirichlet BCs are accounted for in the assembly of the
1020 * RHS; in fact, nothing to be done in this program since we
1021 * prescribe homogeneous BCs).
1022 *
1023
1024 *
1025 *
1026 * @code
1027 *   const typename DoFHandler<dim>::active_cell_iterator
1028 *   neighbor_cell = cell->neighbor(face_no);
1029 *   neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
1030 *  
1031 *   stiffness_matrix_cn = 0;
1032 *   stiffness_matrix_nc = 0;
1033 *   stiffness_matrix_nn = 0;
1034 *   for (unsigned int q = 0; q < n_q_points; ++q)
1035 *   {
1036 *   const double dx = fe_values.JxW(q);
1037 *  
1038 *   for (unsigned int i = 0; i < n_dofs; ++i)
1039 *   {
1040 *   for (unsigned int j = 0; j < n_dofs; ++j)
1041 *   {
1042 *   const Tensor<2, dim> &H_i = discrete_hessians[i][q];
1043 *   const Tensor<2, dim> &H_j = discrete_hessians[j][q];
1044 *  
1045 *   const Tensor<2, dim> &H_i_neigh =
1046 *   discrete_hessians_neigh[face_no][i][q];
1047 *   const Tensor<2, dim> &H_j_neigh =
1048 *   discrete_hessians_neigh[face_no][j][q];
1049 *  
1050 *   stiffness_matrix_cn(i, j) +=
1051 *   scalar_product(H_j_neigh, H_i) * dx;
1052 *   stiffness_matrix_nc(i, j) +=
1053 *   scalar_product(H_j, H_i_neigh) * dx;
1054 *   stiffness_matrix_nn(i, j) +=
1055 *   scalar_product(H_j_neigh, H_i_neigh) * dx;
1056 *   }
1057 *   }
1058 *   }
1059 *  
1060 *   for (unsigned int i = 0; i < n_dofs; ++i)
1061 *   {
1062 *   for (unsigned int j = 0; j < n_dofs; ++j)
1063 *   {
1064 *   matrix(local_dof_indices[i],
1065 *   local_dof_indices_neighbor[j]) +=
1066 *   stiffness_matrix_cn(i, j);
1067 *   matrix(local_dof_indices_neighbor[i],
1068 *   local_dof_indices[j]) +=
1069 *   stiffness_matrix_nc(i, j);
1070 *   matrix(local_dof_indices_neighbor[i],
1071 *   local_dof_indices_neighbor[j]) +=
1072 *   stiffness_matrix_nn(i, j);
1073 *   }
1074 *   }
1075 *  
1076 *   } // boundary check
1077 *   } // for face
1078 *  
1079 * @endcode
1080 *
1081 * We now compute and add the interactions of the degrees of freedom of
1082 * a neighboring cells with those of another neighboring cell (this is
1083 * where we need the extended sparsity pattern).
1084 *
1085 * @code
1086 *   for (unsigned int face_no = 0; face_no < cell->n_faces() - 1; ++face_no)
1087 *   {
1088 *   const typename DoFHandler<dim>::face_iterator face =
1089 *   cell->face(face_no);
1090 *  
1091 *   const bool at_boundary = face->at_boundary();
1092 *   if (!at_boundary)
1093 *   { // nothing to be done if boundary face (the liftings of the
1094 * @endcode
1095 *
1096 * Dirichlet BCs are accounted for in the assembly of the RHS;
1097 * in fact, nothing to be done in this program since we
1098 * prescribe homogeneous BCs)
1099 *
1100
1101 *
1102 *
1103
1104 *
1105 *
1106 * @code
1107 *   for (unsigned int face_no_2 = face_no + 1;
1108 *   face_no_2 < cell->n_faces();
1109 *   ++face_no_2)
1110 *   {
1111 *   const typename DoFHandler<dim>::face_iterator face_2 =
1112 *   cell->face(face_no_2);
1113 *  
1114 *   const bool at_boundary_2 = face_2->at_boundary();
1115 *   if (!at_boundary_2)
1116 *   {
1117 *   const typename DoFHandler<dim>::active_cell_iterator
1118 *   neighbor_cell = cell->neighbor(face_no);
1119 *   neighbor_cell->get_dof_indices(
1120 *   local_dof_indices_neighbor);
1121 *   const typename DoFHandler<dim>::active_cell_iterator
1122 *   neighbor_cell_2 = cell->neighbor(face_no_2);
1123 *   neighbor_cell_2->get_dof_indices(
1124 *   local_dof_indices_neighbor_2);
1125 *  
1126 *   stiffness_matrix_n1n2 = 0;
1127 *   stiffness_matrix_n2n1 = 0;
1128 *  
1129 *   for (unsigned int q = 0; q < n_q_points; ++q)
1130 *   {
1131 *   const double dx = fe_values.JxW(q);
1132 *  
1133 *   for (unsigned int i = 0; i < n_dofs; ++i)
1134 *   for (unsigned int j = 0; j < n_dofs; ++j)
1135 *   {
1136 *   const Tensor<2, dim> &H_i_neigh =
1137 *   discrete_hessians_neigh[face_no][i][q];
1138 *   const Tensor<2, dim> &H_j_neigh =
1139 *   discrete_hessians_neigh[face_no][j][q];
1140 *  
1141 *   const Tensor<2, dim> &H_i_neigh2 =
1142 *   discrete_hessians_neigh[face_no_2][i][q];
1143 *   const Tensor<2, dim> &H_j_neigh2 =
1144 *   discrete_hessians_neigh[face_no_2][j][q];
1145 *  
1146 *   stiffness_matrix_n1n2(i, j) +=
1147 *   scalar_product(H_j_neigh2, H_i_neigh) * dx;
1148 *   stiffness_matrix_n2n1(i, j) +=
1149 *   scalar_product(H_j_neigh, H_i_neigh2) * dx;
1150 *   }
1151 *   }
1152 *  
1153 *   for (unsigned int i = 0; i < n_dofs; ++i)
1154 *   for (unsigned int j = 0; j < n_dofs; ++j)
1155 *   {
1156 *   matrix(local_dof_indices_neighbor[i],
1157 *   local_dof_indices_neighbor_2[j]) +=
1158 *   stiffness_matrix_n1n2(i, j);
1159 *   matrix(local_dof_indices_neighbor_2[i],
1160 *   local_dof_indices_neighbor[j]) +=
1161 *   stiffness_matrix_n2n1(i, j);
1162 *   }
1163 *   } // boundary check face_2
1164 *   } // for face_2
1165 *   } // boundary check face_1
1166 *   } // for face_1
1167 *  
1168 *  
1169 * @endcode
1170 *
1171 * Finally, we compute and add the two penalty terms.
1172 *
1173 * @code
1174 *   for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1175 *   {
1176 *   const typename DoFHandler<dim>::face_iterator face =
1177 *   cell->face(face_no);
1178 *  
1179 *   const double mesh_inv = 1.0 / face->diameter(); // h_e^{-1}
1180 *   const double mesh3_inv =
1181 *   1.0 / Utilities::fixed_power<3>(face->diameter()); // h_e^{-3}
1182 *  
1183 *   fe_face.reinit(cell, face_no);
1184 *  
1185 *   ip_matrix_cc = 0; // filled in any case (boundary or interior face)
1186 *  
1187 *   const bool at_boundary = face->at_boundary();
1188 *   if (at_boundary)
1189 *   {
1190 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1191 *   {
1192 *   const double dx = fe_face.JxW(q);
1193 *  
1194 *   for (unsigned int i = 0; i < n_dofs; ++i)
1195 *   for (unsigned int j = 0; j < n_dofs; ++j)
1196 *   {
1197 *   ip_matrix_cc(i, j) += penalty_jump_grad * mesh_inv *
1198 *   fe_face.shape_grad(j, q) *
1199 *   fe_face.shape_grad(i, q) * dx;
1200 *   ip_matrix_cc(i, j) += penalty_jump_val * mesh3_inv *
1201 *   fe_face.shape_value(j, q) *
1202 *   fe_face.shape_value(i, q) * dx;
1203 *   }
1204 *   }
1205 *   }
1206 *   else
1207 *   { // interior face
1208 *  
1209 *   const typename DoFHandler<dim>::active_cell_iterator
1210 *   neighbor_cell = cell->neighbor(face_no);
1211 *   const unsigned int face_no_neighbor =
1212 *   cell->neighbor_of_neighbor(face_no);
1213 *  
1214 * @endcode
1215 *
1216 * In the next step, we need to have a global way to compare the
1217 * cells in order to not calculate the same jump term twice:
1218 *
1219 * @code
1220 *   if (neighbor_cell->id() < cell->id())
1221 *   continue; // skip this face (already considered)
1222 *   else
1223 *   {
1224 *   fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1225 *   neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
1226 *  
1227 *   ip_matrix_cn = 0;
1228 *   ip_matrix_nc = 0;
1229 *   ip_matrix_nn = 0;
1230 *  
1231 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1232 *   {
1233 *   const double dx = fe_face.JxW(q);
1234 *  
1235 *   for (unsigned int i = 0; i < n_dofs; ++i)
1236 *   {
1237 *   for (unsigned int j = 0; j < n_dofs; ++j)
1238 *   {
1239 *   ip_matrix_cc(i, j) +=
1240 *   penalty_jump_grad * mesh_inv *
1241 *   fe_face.shape_grad(j, q) *
1242 *   fe_face.shape_grad(i, q) * dx;
1243 *   ip_matrix_cc(i, j) +=
1244 *   penalty_jump_val * mesh3_inv *
1245 *   fe_face.shape_value(j, q) *
1246 *   fe_face.shape_value(i, q) * dx;
1247 *  
1248 *   ip_matrix_cn(i, j) -=
1249 *   penalty_jump_grad * mesh_inv *
1250 *   fe_face_neighbor.shape_grad(j, q) *
1251 *   fe_face.shape_grad(i, q) * dx;
1252 *   ip_matrix_cn(i, j) -=
1253 *   penalty_jump_val * mesh3_inv *
1254 *   fe_face_neighbor.shape_value(j, q) *
1255 *   fe_face.shape_value(i, q) * dx;
1256 *  
1257 *   ip_matrix_nc(i, j) -=
1258 *   penalty_jump_grad * mesh_inv *
1259 *   fe_face.shape_grad(j, q) *
1260 *   fe_face_neighbor.shape_grad(i, q) * dx;
1261 *   ip_matrix_nc(i, j) -=
1262 *   penalty_jump_val * mesh3_inv *
1263 *   fe_face.shape_value(j, q) *
1264 *   fe_face_neighbor.shape_value(i, q) * dx;
1265 *  
1266 *   ip_matrix_nn(i, j) +=
1267 *   penalty_jump_grad * mesh_inv *
1268 *   fe_face_neighbor.shape_grad(j, q) *
1269 *   fe_face_neighbor.shape_grad(i, q) * dx;
1270 *   ip_matrix_nn(i, j) +=
1271 *   penalty_jump_val * mesh3_inv *
1272 *   fe_face_neighbor.shape_value(j, q) *
1273 *   fe_face_neighbor.shape_value(i, q) * dx;
1274 *   }
1275 *   }
1276 *   }
1277 *   } // face not visited yet
1278 *  
1279 *   } // boundary check
1280 *  
1281 *   for (unsigned int i = 0; i < n_dofs; ++i)
1282 *   {
1283 *   for (unsigned int j = 0; j < n_dofs; ++j)
1284 *   {
1285 *   matrix(local_dof_indices[i], local_dof_indices[j]) +=
1286 *   ip_matrix_cc(i, j);
1287 *   }
1288 *   }
1289 *  
1290 *   if (!at_boundary)
1291 *   {
1292 *   for (unsigned int i = 0; i < n_dofs; ++i)
1293 *   {
1294 *   for (unsigned int j = 0; j < n_dofs; ++j)
1295 *   {
1296 *   matrix(local_dof_indices[i],
1297 *   local_dof_indices_neighbor[j]) +=
1298 *   ip_matrix_cn(i, j);
1299 *   matrix(local_dof_indices_neighbor[i],
1300 *   local_dof_indices[j]) += ip_matrix_nc(i, j);
1301 *   matrix(local_dof_indices_neighbor[i],
1302 *   local_dof_indices_neighbor[j]) +=
1303 *   ip_matrix_nn(i, j);
1304 *   }
1305 *   }
1306 *   }
1307 *  
1308 *   } // for face
1309 *   } // for cell
1310 *   }
1311 *  
1312 *  
1313 *  
1314 * @endcode
1315 *
1316 *
1317 * <a name="step_82-BiLaplacianLDGLiftassemble_rhs"></a>
1318 * <h4>BiLaplacianLDGLift::assemble_rhs</h4>
1319 *
1320
1321 *
1322 * This function assemble the right-hand side of the system. Since we consider
1323 * homogeneous Dirichlet boundary conditions, imposed weakly in the bilinear
1324 * form using the Nitsche approach, it only involves the contribution of the
1325 * forcing term @f$\int_{\Omega}fv_h@f$.
1326 *
1327 * @code
1328 *   template <int dim>
1329 *   void BiLaplacianLDGLift<dim>::assemble_rhs()
1330 *   {
1331 *   rhs = 0;
1332 *  
1333 *   const QGauss<dim> quad(fe.degree + 1);
1334 *   FEValues<dim> fe_values(
1336 *  
1337 *   const unsigned int n_dofs = fe_values.dofs_per_cell;
1338 *   const unsigned int n_quad_pts = quad.size();
1339 *  
1340 *   const RightHandSide<dim> right_hand_side;
1341 *  
1342 *   Vector<double> local_rhs(n_dofs);
1343 *   std::vector<types::global_dof_index> local_dof_indices(n_dofs);
1344 *  
1345 *   for (const auto &cell : dof_handler.active_cell_iterators())
1346 *   {
1347 *   fe_values.reinit(cell);
1348 *   cell->get_dof_indices(local_dof_indices);
1349 *  
1350 *   local_rhs = 0;
1351 *   for (unsigned int q = 0; q < n_quad_pts; ++q)
1352 *   {
1353 *   const double dx = fe_values.JxW(q);
1354 *  
1355 *   for (unsigned int i = 0; i < n_dofs; ++i)
1356 *   {
1357 *   local_rhs(i) +=
1358 *   right_hand_side.value(fe_values.quadrature_point(q)) *
1359 *   fe_values.shape_value(i, q) * dx;
1360 *   }
1361 *   }
1362 *  
1363 *   for (unsigned int i = 0; i < n_dofs; ++i)
1364 *   rhs(local_dof_indices[i]) += local_rhs(i);
1365 *   }
1366 *   }
1367 *  
1368 *  
1369 *  
1370 * @endcode
1371 *
1372 *
1373 * <a name="step_82-BiLaplacianLDGLiftsolve"></a>
1374 * <h4>BiLaplacianLDGLift::solve</h4>
1375 *
1376
1377 *
1378 * To solve the linear system @f$A\boldsymbol{U}=\boldsymbol{F}@f$,
1379 * we proceed as in @ref step_74 "step-74" and use a direct method. We could
1380 * as well use an iterative method, for instance the conjugate
1381 * gradient method as in @ref step_3 "step-3".
1382 *
1383 * @code
1384 *   template <int dim>
1385 *   void BiLaplacianLDGLift<dim>::solve()
1386 *   {
1387 *   SparseDirectUMFPACK A_direct;
1388 *   A_direct.initialize(matrix);
1389 *   A_direct.vmult(solution, rhs);
1390 *   }
1391 *  
1392 *  
1393 *  
1394 * @endcode
1395 *
1396 *
1397 * <a name="step_82-BiLaplacianLDGLiftcompute_errors"></a>
1398 * <h4>BiLaplacianLDGLift::compute_errors</h4>
1399 *
1400
1401 *
1402 * This function computes the discrete @f$H^2@f$, @f$H^1@f$ and @f$L^2@f$ norms of
1403 * the error @f$u-u_h@f$, where @f$u@f$ is the exact solution and @f$u_h@f$ is
1404 * the approximate solution. See the introduction for the definition
1405 * of the norms.
1406 *
1407 * @code
1408 *   template <int dim>
1409 *   void BiLaplacianLDGLift<dim>::compute_errors()
1410 *   {
1411 *   double error_H2 = 0;
1412 *   double error_H1 = 0;
1413 *   double error_L2 = 0;
1414 *  
1415 *   const QGauss<dim> quad(fe.degree + 1);
1416 *   const QGauss<dim - 1> quad_face(fe.degree + 1);
1417 *  
1418 *   FEValues<dim> fe_values(fe,
1419 *   quad,
1422 *  
1423 *   FEFaceValues<dim> fe_face(fe,
1424 *   quad_face,
1427 *  
1428 *   FEFaceValues<dim> fe_face_neighbor(fe,
1429 *   quad_face,
1431 *  
1432 *   const unsigned int n_q_points = quad.size();
1433 *   const unsigned int n_q_points_face = quad_face.size();
1434 *  
1435 * @endcode
1436 *
1437 * We introduce some variables for the exact solution...
1438 *
1439 * @code
1440 *   const ExactSolution<dim> u_exact;
1441 *  
1442 * @endcode
1443 *
1444 * ...and for the approximate solution:
1445 *
1446 * @code
1447 *   std::vector<double> solution_values_cell(n_q_points);
1448 *   std::vector<Tensor<1, dim>> solution_gradients_cell(n_q_points);
1449 *   std::vector<Tensor<2, dim>> solution_hessians_cell(n_q_points);
1450 *  
1451 *   std::vector<double> solution_values(n_q_points_face);
1452 *   std::vector<double> solution_values_neigh(n_q_points_face);
1453 *   std::vector<Tensor<1, dim>> solution_gradients(n_q_points_face);
1454 *   std::vector<Tensor<1, dim>> solution_gradients_neigh(n_q_points_face);
1455 *  
1456 *   for (const auto &cell : dof_handler.active_cell_iterators())
1457 *   {
1458 *   fe_values.reinit(cell);
1459 *  
1460 *   fe_values.get_function_values(solution, solution_values_cell);
1461 *   fe_values.get_function_gradients(solution, solution_gradients_cell);
1462 *   fe_values.get_function_hessians(solution, solution_hessians_cell);
1463 *  
1464 * @endcode
1465 *
1466 * We first add the <i>bulk</i> terms.
1467 *
1468 * @code
1469 *   for (unsigned int q = 0; q < n_q_points; ++q)
1470 *   {
1471 *   const double dx = fe_values.JxW(q);
1472 *  
1473 *   error_H2 += (u_exact.hessian(fe_values.quadrature_point(q)) -
1474 *   solution_hessians_cell[q])
1475 *   .norm_square() *
1476 *   dx;
1477 *   error_H1 += (u_exact.gradient(fe_values.quadrature_point(q)) -
1478 *   solution_gradients_cell[q])
1479 *   .norm_square() *
1480 *   dx;
1481 *   error_L2 += Utilities::fixed_power<2>(
1482 *   u_exact.value(fe_values.quadrature_point(q)) -
1483 *   solution_values_cell[q]) *
1484 *   dx;
1485 *   } // for quadrature points
1486 *  
1487 * @endcode
1488 *
1489 * We then add the face contributions.
1490 *
1491 * @code
1492 *   for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1493 *   {
1494 *   const typename DoFHandler<dim>::face_iterator face =
1495 *   cell->face(face_no);
1496 *  
1497 *   const double mesh_inv = 1.0 / face->diameter(); // h^{-1}
1498 *   const double mesh3_inv =
1499 *   1.0 / Utilities::fixed_power<3>(face->diameter()); // h^{-3}
1500 *  
1501 *   fe_face.reinit(cell, face_no);
1502 *  
1503 *   fe_face.get_function_values(solution, solution_values);
1504 *   fe_face.get_function_gradients(solution, solution_gradients);
1505 *  
1506 *   const bool at_boundary = face->at_boundary();
1507 *   if (at_boundary)
1508 *   {
1509 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1510 *   {
1511 *   const double dx = fe_face.JxW(q);
1512 *   const double u_exact_q =
1513 *   u_exact.value(fe_face.quadrature_point(q));
1514 *   const Tensor<1, dim> u_exact_grad_q =
1515 *   u_exact.gradient(fe_face.quadrature_point(q));
1516 *  
1517 *   error_H2 +=
1518 *   mesh_inv *
1519 *   (u_exact_grad_q - solution_gradients[q]).norm_square() *
1520 *   dx;
1521 *   error_H2 += mesh3_inv *
1522 *   Utilities::fixed_power<2>(u_exact_q -
1523 *   solution_values[q]) *
1524 *   dx;
1525 *   error_H1 += mesh_inv *
1526 *   Utilities::fixed_power<2>(u_exact_q -
1527 *   solution_values[q]) *
1528 *   dx;
1529 *   }
1530 *   }
1531 *   else
1532 *   { // interior face
1533 *  
1534 *   const typename DoFHandler<dim>::active_cell_iterator
1535 *   neighbor_cell = cell->neighbor(face_no);
1536 *   const unsigned int face_no_neighbor =
1537 *   cell->neighbor_of_neighbor(face_no);
1538 *  
1539 * @endcode
1540 *
1541 * In the next step, we need to have a global way to compare the
1542 * cells in order to not calculate the same jump term twice:
1543 *
1544 * @code
1545 *   if (neighbor_cell->id() < cell->id())
1546 *   continue; // skip this face (already considered)
1547 *   else
1548 *   {
1549 *   fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1550 *  
1551 *   fe_face.get_function_values(solution, solution_values);
1552 *   fe_face_neighbor.get_function_values(solution,
1553 *   solution_values_neigh);
1554 *   fe_face.get_function_gradients(solution,
1555 *   solution_gradients);
1556 *   fe_face_neighbor.get_function_gradients(
1557 *   solution, solution_gradients_neigh);
1558 *  
1559 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1560 *   {
1561 *   const double dx = fe_face.JxW(q);
1562 *  
1563 * @endcode
1564 *
1565 * To compute the jump term, we use the fact that
1566 * @f$\jump{u}=0@f$ and
1567 * @f$\jump{\nabla u}=\mathbf{0}@f$ since @f$u\in
1568 * H^2(\Omega)@f$.
1569 *
1570 * @code
1571 *   error_H2 +=
1572 *   mesh_inv *
1573 *   (solution_gradients_neigh[q] - solution_gradients[q])
1574 *   .norm_square() *
1575 *   dx;
1576 *   error_H2 +=
1577 *   mesh3_inv *
1578 *   Utilities::fixed_power<2>(solution_values_neigh[q] -
1579 *   solution_values[q]) *
1580 *   dx;
1581 *   error_H1 +=
1582 *   mesh_inv *
1583 *   Utilities::fixed_power<2>(solution_values_neigh[q] -
1584 *   solution_values[q]) *
1585 *   dx;
1586 *   }
1587 *   } // face not visited yet
1588 *  
1589 *   } // boundary check
1590 *  
1591 *   } // for face
1592 *  
1593 *   } // for cell
1594 *  
1595 *   error_H2 = std::sqrt(error_H2);
1596 *   error_H1 = std::sqrt(error_H1);
1597 *   error_L2 = std::sqrt(error_L2);
1598 *  
1599 *   std::cout << "DG H2 norm of the error: " << error_H2 << std::endl;
1600 *   std::cout << "DG H1 norm of the error: " << error_H1 << std::endl;
1601 *   std::cout << " L2 norm of the error: " << error_L2 << std::endl;
1602 *   }
1603 *  
1604 *  
1605 *  
1606 * @endcode
1607 *
1608 *
1609 * <a name="step_82-BiLaplacianLDGLiftoutput_results"></a>
1610 * <h4>BiLaplacianLDGLift::output_results</h4>
1611 *
1612
1613 *
1614 * This function, which writes the solution to a vtk file,
1615 * is copied from @ref step_3 "step-3".
1616 *
1617 * @code
1618 *   template <int dim>
1619 *   void BiLaplacianLDGLift<dim>::output_results() const
1620 *   {
1621 *   DataOut<dim> data_out;
1622 *   data_out.attach_dof_handler(dof_handler);
1623 *   data_out.add_data_vector(solution, "solution");
1624 *   data_out.build_patches();
1625 *  
1626 *   std::ofstream output("solution.vtk");
1627 *   data_out.write_vtk(output);
1628 *   }
1629 *  
1630 *  
1631 *  
1632 * @endcode
1633 *
1634 *
1635 * <a name="step_82-BiLaplacianLDGLiftassemble_local_matrix"></a>
1636 * <h4>BiLaplacianLDGLift::assemble_local_matrix</h4>
1637 *
1638
1639 *
1640 * As already mentioned above, this function is used to assemble
1641 * the (local) mass matrices needed for the computations of the
1642 * lifting terms. We reiterate that only the basis functions with
1643 * support on the current cell are considered.
1644 *
1645 * @code
1646 *   template <int dim>
1647 *   void BiLaplacianLDGLift<dim>::assemble_local_matrix(
1648 *   const FEValues<dim> &fe_values_lift,
1649 *   const unsigned int n_q_points,
1650 *   FullMatrix<double> &local_matrix)
1651 *   {
1652 *   const FEValuesExtractors::Tensor<2> tau_ext(0);
1653 *  
1654 *   const unsigned int n_dofs = fe_values_lift.dofs_per_cell;
1655 *  
1656 *   local_matrix = 0;
1657 *   for (unsigned int q = 0; q < n_q_points; ++q)
1658 *   {
1659 *   const double dx = fe_values_lift.JxW(q);
1660 *  
1661 *   for (unsigned int m = 0; m < n_dofs; ++m)
1662 *   for (unsigned int n = 0; n < n_dofs; ++n)
1663 *   {
1664 *   local_matrix(m, n) +=
1665 *   scalar_product(fe_values_lift[tau_ext].value(n, q),
1666 *   fe_values_lift[tau_ext].value(m, q)) *
1667 *   dx;
1668 *   }
1669 *   }
1670 *   }
1671 *  
1672 *  
1673 *  
1674 * @endcode
1675 *
1676 *
1677 * <a name="step_82-BiLaplacianLDGLiftcompute_discrete_hessians"></a>
1678 * <h4>BiLaplacianLDGLift::compute_discrete_hessians</h4>
1679 *
1680
1681 *
1682 * This function is the main novelty of this program. It computes
1683 * the discrete Hessian @f$H_h(\varphi)@f$ for all the basis functions
1684 * @f$\varphi@f$ of @f$\mathbb{V}_h@f$ supported on the current cell and
1685 * those supported on a neighboring cell. The first argument
1686 * indicates the current cell (referring to the global DoFHandler
1687 * object), while the other two arguments are output variables that
1688 * are filled by this function.
1689 *
1690
1691 *
1692 * In the following, we need to evaluate finite element shape
1693 * functions for the `fe_lift` finite element on the current
1694 * cell. Like for example in @ref step_61 "step-61", this "lift" space is defined
1695 * on every cell individually; as a consequence, there is no global
1696 * DoFHandler associated with this because we simply have no need
1697 * for such a DoFHandler. That leaves the question of what we should
1698 * initialize the FEValues and FEFaceValues objects with when we ask
1699 * them to evaluate shape functions of `fe_lift` on a concrete
1700 * cell. If we simply provide the first argument to this function,
1701 * `cell`, to FEValues::reinit(), we will receive an error message
1702 * that the given `cell` belongs to a DoFHandler that has a
1703 * different finite element associated with it than the `fe_lift`
1704 * object we want to evaluate. Fortunately, there is a relatively
1705 * easy solution: We can call FEValues::reinit() with a cell that
1706 * points into a triangulation -- the same cell, but not associated
1707 * with a DoFHandler, and consequently no finite element space. In
1708 * that case, FEValues::reinit() will skip the check that would
1709 * otherwise lead to an error message. All we have to do is to convert
1710 * the DoFHandler cell iterator into a Triangulation cell iterator;
1711 * see the first couple of lines of the function below to see how
1712 * this can be done.
1713 *
1714 * @code
1715 *   template <int dim>
1716 *   void BiLaplacianLDGLift<dim>::compute_discrete_hessians(
1717 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1718 *   std::vector<std::vector<Tensor<2, dim>>> &discrete_hessians,
1719 *   std::vector<std::vector<std::vector<Tensor<2, dim>>>>
1720 *   &discrete_hessians_neigh)
1721 *   {
1722 *   const typename Triangulation<dim>::cell_iterator cell_lift =
1723 *   static_cast<typename Triangulation<dim>::cell_iterator>(cell);
1724 *  
1725 *   const QGauss<dim> quad(fe.degree + 1);
1726 *   const QGauss<dim - 1> quad_face(fe.degree + 1);
1727 *  
1728 *   const unsigned int n_q_points = quad.size();
1729 *   const unsigned int n_q_points_face = quad_face.size();
1730 *  
1731 * @endcode
1732 *
1733 * The information we need from the basis functions of
1734 * @f$\mathbb{V}_h@f$: <code>fe_values</code> is needed to add
1735 * the broken Hessian part of the discrete Hessian, while
1736 * <code>fe_face</code> and <code>fe_face_neighbor</code>
1737 * are used to compute the right-hand sides for the local
1738 * problems.
1739 *
1740 * @code
1741 *   FEValues<dim> fe_values(fe, quad, update_hessians | update_JxW_values);
1742 *  
1743 *   FEFaceValues<dim> fe_face(
1745 *  
1746 *   FEFaceValues<dim> fe_face_neighbor(
1748 *  
1749 *   const unsigned int n_dofs = fe_values.dofs_per_cell;
1750 *  
1751 * @endcode
1752 *
1753 * The information needed from the basis functions
1754 * of the finite element space for the lifting terms:
1755 * <code>fe_values_lift</code> is used for the (local)
1756 * mass matrix (see @f$\boldsymbol{M}_c@f$ in the introduction),
1757 * while <code>fe_face_lift</code> is used to compute the
1758 * right-hand sides (see @f$\boldsymbol{G}_c@f$ for @f$b_e@f$).
1759 *
1760 * @code
1761 *   FEValues<dim> fe_values_lift(fe_lift,
1762 *   quad,
1764 *  
1765 *   FEFaceValues<dim> fe_face_lift(
1766 *   fe_lift, quad_face, update_values | update_gradients | update_JxW_values);
1767 *  
1768 *   const FEValuesExtractors::Tensor<2> tau_ext(0);
1769 *  
1770 *   const unsigned int n_dofs_lift = fe_values_lift.dofs_per_cell;
1771 *   FullMatrix<double> local_matrix_lift(n_dofs_lift, n_dofs_lift);
1772 *  
1773 *   Vector<double> local_rhs_re(n_dofs_lift), local_rhs_be(n_dofs_lift),
1774 *   coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift), coeffs_tmp(n_dofs_lift);
1775 *  
1776 *   SolverControl solver_control(1000, 1e-12);
1777 *   SolverCG<Vector<double>> solver(solver_control);
1778 *  
1779 *   double factor_avg; // 0.5 for interior faces, 1.0 for boundary faces
1780 *  
1781 *   fe_values.reinit(cell);
1782 *   fe_values_lift.reinit(cell_lift);
1783 *  
1784 * @endcode
1785 *
1786 * We start by assembling the (local) mass matrix used for the computation
1787 * of the lifting terms @f$r_e@f$ and @f$b_e@f$.
1788 *
1789 * @code
1790 *   assemble_local_matrix(fe_values_lift, n_q_points, local_matrix_lift);
1791 *  
1792 *   for (unsigned int i = 0; i < n_dofs; ++i)
1793 *   for (unsigned int q = 0; q < n_q_points; ++q)
1794 *   {
1795 *   discrete_hessians[i][q] = 0;
1796 *  
1797 *   for (unsigned int face_no = 0;
1798 *   face_no < discrete_hessians_neigh.size();
1799 *   ++face_no)
1800 *   {
1801 *   discrete_hessians_neigh[face_no][i][q] = 0;
1802 *   }
1803 *   }
1804 *  
1805 * @endcode
1806 *
1807 * In this loop, we compute the discrete Hessian at each quadrature point
1808 * @f$x_q@f$ of <code>cell</code> for each basis function supported on
1809 * <code>cell</code>, namely we fill-in the variable
1810 * <code>discrete_hessians[i][q]</code>. For the lifting terms, we need to
1811 * add the contribution of all the faces of <code>cell</code>.
1812 *
1813 * @code
1814 *   for (unsigned int i = 0; i < n_dofs; ++i)
1815 *   {
1816 *   coeffs_re = 0;
1817 *   coeffs_be = 0;
1818 *  
1819 *   for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1820 *   {
1821 *   const typename DoFHandler<dim>::face_iterator face =
1822 *   cell->face(face_no);
1823 *  
1824 *   const bool at_boundary = face->at_boundary();
1825 *  
1826 * @endcode
1827 *
1828 * Recall that by convention, the average of a function across a
1829 * boundary face @f$e@f$ reduces to the trace of the function on the
1830 * only element adjacent to @f$e@f$, namely there is no factor
1831 * @f$\frac{1}{2}@f$. We distinguish between the two cases (the current
1832 * face lies in the interior or on the boundary of the domain) using
1833 * the variable <code>factor_avg</code>.
1834 *
1835 * @code
1836 *   factor_avg = 0.5;
1837 *   if (at_boundary)
1838 *   {
1839 *   factor_avg = 1.0;
1840 *   }
1841 *  
1842 *   fe_face.reinit(cell, face_no);
1843 *   fe_face_lift.reinit(cell_lift, face_no);
1844 *  
1845 *   local_rhs_re = 0;
1846 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1847 *   {
1848 *   const double dx = fe_face_lift.JxW(q);
1849 *   const Tensor<1, dim> normal = fe_face.normal_vector(
1850 *   q); // same as fe_face_lift.normal_vector(q)
1851 *  
1852 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1853 *   {
1854 *   local_rhs_re(m) +=
1855 *   factor_avg *
1856 *   (fe_face_lift[tau_ext].value(m, q) * normal) *
1857 *   fe_face.shape_grad(i, q) * dx;
1858 *   }
1859 *   }
1860 *  
1861 * @endcode
1862 *
1863 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
1864 * introduced in the comments about the implementation of the
1865 * lifting @f$b_e@f$ in the case
1866 * @f$\varphi=\varphi^c@f$.
1867 *
1868 * @code
1869 *   local_rhs_be = 0;
1870 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1871 *   {
1872 *   const double dx = fe_face_lift.JxW(q);
1873 *   const Tensor<1, dim> normal = fe_face.normal_vector(
1874 *   q); // same as fe_face_lift.normal_vector(q)
1875 *  
1876 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1877 *   {
1878 *   local_rhs_be(m) += factor_avg *
1879 *   fe_face_lift[tau_ext].divergence(m, q) *
1880 *   normal * fe_face.shape_value(i, q) * dx;
1881 *   }
1882 *   }
1883 *  
1884 *   coeffs_tmp = 0;
1885 *   solver.solve(local_matrix_lift,
1886 *   coeffs_tmp,
1887 *   local_rhs_re,
1888 *   PreconditionIdentity());
1889 *   coeffs_re += coeffs_tmp;
1890 *  
1891 *   coeffs_tmp = 0;
1892 *   solver.solve(local_matrix_lift,
1893 *   coeffs_tmp,
1894 *   local_rhs_be,
1895 *   PreconditionIdentity());
1896 *   coeffs_be += coeffs_tmp;
1897 *  
1898 *   } // for face
1899 *  
1900 *   for (unsigned int q = 0; q < n_q_points; ++q)
1901 *   {
1902 *   discrete_hessians[i][q] += fe_values.shape_hessian(i, q);
1903 *  
1904 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1905 *   {
1906 *   discrete_hessians[i][q] -=
1907 *   coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
1908 *   }
1909 *  
1910 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1911 *   {
1912 *   discrete_hessians[i][q] +=
1913 *   coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
1914 *   }
1915 *   }
1916 *   } // for dof i
1917 *  
1918 *  
1919 *  
1920 * @endcode
1921 *
1922 * In this loop, we compute the discrete Hessian at each quadrature point
1923 * @f$x_q@f$ of <code>cell</code> for each basis function supported on a
1924 * neighboring <code>neighbor_cell</code> of <code>cell</code>, namely we
1925 * fill-in the variable <code>discrete_hessians_neigh[face_no][i][q]</code>.
1926 * For the lifting terms, we only need to add the contribution of the
1927 * face adjacent to <code>cell</code> and <code>neighbor_cell</code>.
1928 *
1929 * @code
1930 *   for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1931 *   {
1932 *   const typename DoFHandler<dim>::face_iterator face =
1933 *   cell->face(face_no);
1934 *  
1935 *   const bool at_boundary = face->at_boundary();
1936 *  
1937 *   if (!at_boundary)
1938 *   {
1939 * @endcode
1940 *
1941 * For non-homogeneous Dirichlet BCs, we would need to
1942 * compute the lifting of the prescribed BC (see the
1943 * "Possible Extensions" section for more details).
1944 *
1945
1946 *
1947 *
1948 * @code
1950 *   neighbor_cell = cell->neighbor(face_no);
1951 *   const unsigned int face_no_neighbor =
1952 *   cell->neighbor_of_neighbor(face_no);
1953 *   fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1954 *  
1955 *   for (unsigned int i = 0; i < n_dofs; ++i)
1956 *   {
1957 *   coeffs_re = 0;
1958 *   coeffs_be = 0;
1959 *  
1960 *   fe_face_lift.reinit(cell_lift, face_no);
1961 *  
1962 *   local_rhs_re = 0;
1963 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1964 *   {
1965 *   const double dx = fe_face_lift.JxW(q);
1966 *   const Tensor<1, dim> normal =
1967 *   fe_face_neighbor.normal_vector(q);
1968 *  
1969 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1970 *   {
1971 *   local_rhs_re(m) +=
1972 *   0.5 * (fe_face_lift[tau_ext].value(m, q) * normal) *
1973 *   fe_face_neighbor.shape_grad(i, q) * dx;
1974 *   }
1975 *   }
1976 *  
1977 * @endcode
1978 *
1979 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
1980 * introduced in the comments about the implementation of the
1981 * lifting @f$b_e@f$ in the case
1982 * @f$\varphi=\varphi^n@f$.
1983 *
1984 * @code
1985 *   local_rhs_be = 0;
1986 *   for (unsigned int q = 0; q < n_q_points_face; ++q)
1987 *   {
1988 *   const double dx = fe_face_lift.JxW(q);
1989 *   const Tensor<1, dim> normal =
1990 *   fe_face_neighbor.normal_vector(q);
1991 *  
1992 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
1993 *   {
1994 *   local_rhs_be(m) +=
1995 *   0.5 * fe_face_lift[tau_ext].divergence(m, q) *
1996 *   normal * fe_face_neighbor.shape_value(i, q) * dx;
1997 *   }
1998 *   }
1999 *  
2000 *   solver.solve(local_matrix_lift,
2001 *   coeffs_re,
2002 *   local_rhs_re,
2003 *   PreconditionIdentity());
2004 *   solver.solve(local_matrix_lift,
2005 *   coeffs_be,
2006 *   local_rhs_be,
2007 *   PreconditionIdentity());
2008 *  
2009 *   for (unsigned int q = 0; q < n_q_points; ++q)
2010 *   {
2011 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
2012 *   {
2013 *   discrete_hessians_neigh[face_no][i][q] -=
2014 *   coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
2015 *   }
2016 *  
2017 *   for (unsigned int m = 0; m < n_dofs_lift; ++m)
2018 *   {
2019 *   discrete_hessians_neigh[face_no][i][q] +=
2020 *   coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
2021 *   }
2022 *   }
2023 *  
2024 *   } // for dof i
2025 *   } // boundary check
2026 *   } // for face
2027 *   }
2028 *  
2029 *  
2030 *  
2031 * @endcode
2032 *
2033 *
2034 * <a name="step_82-BiLaplacianLDGLiftrun"></a>
2035 * <h4>BiLaplacianLDGLift::run</h4>
2036 *
2037 * @code
2038 *   template <int dim>
2039 *   void BiLaplacianLDGLift<dim>::run()
2040 *   {
2041 *   make_grid();
2042 *  
2043 *   setup_system();
2044 *   assemble_system();
2045 *  
2046 *   solve();
2047 *  
2048 *   compute_errors();
2049 *   output_results();
2050 *   }
2051 *  
2052 *   } // namespace Step82
2053 *  
2054 *  
2055 *  
2056 * @endcode
2057 *
2058 *
2059 * <a name="step_82-Thecodemaincodefunction"></a>
2060 * <h3>The <code>main</code> function</h3>
2061 *
2062
2063 *
2064 * This is the <code>main</code> function. We define here the number of mesh
2065 * refinements, the polynomial degree for the two finite element spaces
2066 * (for the solution and the two liftings) and the two penalty coefficients.
2067 * We can also change the dimension to run the code in 3d.
2068 *
2069 * @code
2070 *   int main()
2071 *   {
2072 *   try
2073 *   {
2074 *   const unsigned int n_ref = 3; // number of mesh refinements
2075 *  
2076 *   const unsigned int degree =
2077 *   2; // FE degree for u_h and the two lifting terms
2078 *  
2079 *   const double penalty_grad =
2080 *   1.0; // penalty coefficient for the jump of the gradients
2081 *   const double penalty_val =
2082 *   1.0; // penalty coefficient for the jump of the values
2083 *  
2084 *   Step82::BiLaplacianLDGLift<2> problem(n_ref,
2085 *   degree,
2086 *   penalty_grad,
2087 *   penalty_val);
2088 *  
2089 *   problem.run();
2090 *   }
2091 *   catch (std::exception &exc)
2092 *   {
2093 *   std::cerr << std::endl
2094 *   << std::endl
2095 *   << "----------------------------------------------------"
2096 *   << std::endl;
2097 *   std::cerr << "Exception on processing: " << std::endl
2098 *   << exc.what() << std::endl
2099 *   << "Aborting!" << std::endl
2100 *   << "----------------------------------------------------"
2101 *   << std::endl;
2102 *   return 1;
2103 *   }
2104 *   catch (...)
2105 *   {
2106 *   std::cerr << std::endl
2107 *   << std::endl
2108 *   << "----------------------------------------------------"
2109 *   << std::endl;
2110 *   std::cerr << "Unknown exception!" << std::endl
2111 *   << "Aborting!" << std::endl
2112 *   << "----------------------------------------------------"
2113 *   << std::endl;
2114 *   return 1;
2115 *   }
2116 *  
2117 *   return 0;
2118 *   }
2119 * @endcode
2120<a name="step_82-Results"></a><h1>Results</h1>
2121
2122
2123
2124When running the program, the sparsity pattern is written to an svg file, the solution is written to a vtk file, and some results are printed to the console. With the current setup, the output should read
2125
2126@code
2127
2128Number of active cells: 64
2129Number of degrees of freedom: 576
2130Assembling the system.............
2131Done.
2132DG H2 norm of the error: 0.0151063
2133DG H1 norm of the error: 0.000399747
2134 L2 norm of the error: 5.33856e-05
2135
2136@endcode
2137
2138This corresponds to the bi-Laplacian problem with the manufactured solution mentioned above for @f$d=2@f$, 3 refinements of the mesh, degree @f$k=2@f$, and @f$\gamma_0=\gamma_1=1@f$ for the penalty coefficients. By changing the number of refinements, we get the following results:
2139
2140<table align="center" class="doxtable">
2141 <tr>
2142 <th>n_ref</th>
2143 <th>n_cells</th>
2144 <th>n_dofs</th>
2145 <th>error H2 </th>
2146 <th>rate</th>
2147 <th>error H1</th>
2148 <th>rate</th>
2149 <th>error L2</th>
2150 <th>rate</th>
2151 </tr>
2152 <tr>
2153 <td align="center">1</td>
2154 <td align="right">4</td>
2155 <td align="right">36</td>
2156 <td align="center">5.651e-02</td>
2157 <td align="center">--</td>
2158 <td align="center">3.366e-03</td>
2159 <td align="center">--</td>
2160 <td align="center">3.473e-04</td>
2161 <td align="center">--</td>
2162 </tr>
2163 <tr>
2164 <td align="center">2</td>
2165 <td align="right">16</td>
2166 <td align="right">144</td>
2167 <td align="center">3.095e-02</td>
2168 <td align="center">0.87</td>
2169 <td align="center">1.284e-03</td>
2170 <td align="center">1.39</td>
2171 <td align="center">1.369e-04</td>
2172 <td align="center">1.34</td>
2173 </tr>
2174 <tr>
2175 <td align="center">3</td>
2176 <td align="right">64</td>
2177 <td align="right">576</td>
2178 <td align="center">1.511e-02</td>
2179 <td align="center">1.03</td>
2180 <td align="center">3.997e-04</td>
2181 <td align="center">1.68</td>
2182 <td align="center">5.339e-05</td>
2183 <td align="center">1.36</td>
2184 </tr>
2185 <tr>
2186 <td align="center">4</td>
2187 <td align="right">256</td>
2188 <td align="right">2304</td>
2189 <td align="center">7.353e-03</td>
2190 <td align="center">1.04</td>
2191 <td align="center">1.129e-04</td>
2192 <td align="center">1.82</td>
2193 <td align="center">1.691e-05</td>
2194 <td align="center">1.66</td>
2195 </tr>
2196 <tr>
2197 <td align="center">5</td>
2198 <td align="right">1024</td>
2199 <td align="right">9216</td>
2200 <td align="center">3.609e-03</td>
2201 <td align="center">1.03</td>
2202 <td align="center">3.024e-05</td>
2203 <td align="center">1.90</td>
2204 <td align="center">4.789e-06</td>
2205 <td align="center">1.82</td>
2206 </tr>
2207 <tr>
2208 <td align="center">6</td>
2209 <td align="right">4096</td>
2210 <td align="right">36864</td>
2211 <td align="center">1.785e-03</td>
2212 <td align="center">1.02</td>
2213 <td align="center">7.850e-06</td>
2214 <td align="center">1.95</td>
2215 <td align="center">1.277e-06</td>
2216 <td align="center">1.91</td>
2217 </tr>
2218</table>
2219
2220This matches the expected optimal convergence rates for the @f$H^2@f$ and
2221@f$H^1@f$ norms, but is sub-optimal for the @f$L_2@f$ norm. Incidentally, this
2222also matches the results seen in @ref step_47 "step-47" when using polynomial degree
2223@f$k=2@f$.
2224
2225Indeed, just like in @ref step_47 "step-47", we can regain the optimal convergence
2226order if we set the polynomial degree of the finite elements to @f$k=3@f$
2227or higher. Here are the numbers for @f$k=3@f$:
2228
2229<table align="center" class="doxtable">
2230 <tr> <th> n_ref </th> <th> n_cells </th> <th> n_dofs </th> <th> error H2 </th> <th> rate </th> <th> error H1 </th> <th> rate </th> <th> error L2 </th> <th> rate</th> </tr>
2231 <tr> <td> 1 </td> <td> 4 </td> <td> 36 </td> <td> 1.451e-02 </td> <td> -- </td> <td> 5.494e-04 </td> <td> -- </td> <td> 3.035e-05 </td> <td> --</td> </tr>
2232 <tr> <td> 2 </td> <td> 16 </td> <td> 144 </td> <td> 3.565e-03 </td> <td> 2.02 </td> <td> 6.870e-05 </td> <td> 3.00 </td> <td> 2.091e-06 </td> <td> 3.86</td> </tr>
2233 <tr> <td> 3 </td> <td> 64 </td> <td> 576 </td> <td> 8.891e-04 </td> <td> 2.00 </td> <td> 8.584e-06 </td> <td> 3.00 </td> <td> 1.352e-07 </td> <td> 3.95</td> </tr>
2234 <tr> <td> 4 </td> <td> 256 </td> <td> 2304 </td> <td> 2.223e-04 </td> <td> 2.00 </td> <td> 1.073e-06 </td> <td> 3.00 </td> <td> 8.594e-09 </td> <td> 3.98</td> </tr>
2235 <tr> <td> 5 </td> <td> 1024 </td> <td> 9216 </td> <td> 5.560e-05 </td> <td> 2.00 </td> <td> 1.341e-07 </td> <td> 3.00 </td> <td> 5.418e-10 </td> <td> 3.99</td> </tr>
2236 <tr> <td> 6 </td> <td> 4096 </td> <td> 36864 </td> <td> 1.390e-05 </td> <td> 2.00 </td> <td> 1.676e-08 </td> <td> 3.00 </td> <td> 3.245e-11 </td> <td> 4.06</td> </tr>
2237</table>
2238
2239
2240<a name="step_82-Possibleextensions"></a><h3>Possible extensions</h3>
2241
2242
2243The code can be easily adapted to deal with the following cases:
2244
2245<ol>
2246 <li>Non-homogeneous Dirichlet boundary conditions on (part of) the boundary @f$\partial \Omega@f$ of @f$\Omega@f$.</li>
2247 <li>Hanging-nodes (proceed as in @ref step_14 "step-14" to not visit a sub-face twice when computing the lifting terms in <code>compute_discrete_hessians</code> and the penalty terms in <code>assemble_matrix</code>).</li>
2248 <li>LDG method for the Poisson problem (use the discrete gradient consisting of the broken gradient and the lifting of the jump of @f$u_h@f$).</li>
2249</ol>
2250
2251We give below additional details for the first of these points.
2252
2253
2254<a name="step_82-NonhomogeneousDirichletboundaryconditions"></a><h4>Non-homogeneous Dirichlet boundary conditions</h4>
2255
2256If we prescribe non-homogeneous Dirichlet conditions, say
2257@f[
2258\nabla u=\mathbf{g} \quad \mbox{and} \quad u=g \qquad \mbox{on } \partial \Omega,
2259@f]
2260then the right-hand side @f$\boldsymbol{F}@f$ of the linear system needs to be modified as follows
2261@f[
2262F_i:=\int_{\Omega}f\varphi_i-\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}r_e(\mathbf{g}):H_h(\varphi_i)+\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}b_e(g):H_h(\varphi_i)+\gamma_1\sum_{e\in\mathcal{E}_h^b}h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0\sum_{e\in\mathcal{E}_h^b}h_e^{-3}\int_e g\varphi_i, \qquad 1\leq i \leq N_h.
2263@f]
2264Note that for any given index @f$i@f$, many of the terms are zero. Indeed, for @f$e\in \mathcal{E}_h^b@f$ we have @f${\rm supp}\,(r_e(\mathbf{g}))={\rm supp}\,(b_e(g))=K@f$, where @f$K@f$ is the element for which @f$e\subset\partial K@f$. Therefore, the liftings @f$r_e(\mathbf{g})@f$ and @f$b_e(g)@f$ contribute to @f$F_i@f$ only if @f$\varphi_i@f$ has support on @f$K@f$ or a neighbor of @f$K@f$. In other words, when integrating on a cell @f$K@f$, we need to add
2265@f[
2266\int_{K}f\varphi_i+\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)+\gamma_1h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0h_e^{-3}\int_e g\varphi_i\right]
2267@f]
2268to @f$F_i@f$ for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on @f$K@f$ and
2269@f[
2270\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)\right]
2271@f]
2272to @f$F_i@f$ for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on a neighbor of @f$K@f$.
2273
2274@note
2275Note that we can easily consider the case where Dirichlet boundary conditions are imposed only on a subset @f$\emptyset\neq\Gamma_D\subset\partial \Omega@f$. In this case, we simply need to replace @f$\mathcal{E}_h^b@f$ by @f$\mathcal{E}_h^D\subset\mathcal{E}_h^b@f$ consisting of the faces belonging to @f$\Gamma_D@f$. This also affects the matrix @f$A@f$ (simply replace @f$\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^b@f$ by @f$\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^D@f$).
2276 *
2277 *
2278<a name="step_82-PlainProg"></a>
2279<h1> The plain program</h1>
2280@include "step-82.cc"
2281*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
const Quadrature< dim > quadrature
Definition fe_values.h:172
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual SymmetricTensor< 2, dim, RangeNumberType > hessian(const Point< dim > &p, const unsigned int component=0) const
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
void initialize(const SparsityPattern &sparsity_pattern)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
@ 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()
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
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
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
constexpr T fixed_power(const T t)
Definition utilities.h:942
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)
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)