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