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