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