Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
step-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 *
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>
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
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 * {
1105 * neighbor_cell = cell->neighbor(face_no);
1106 * neighbor_cell->get_dof_indices(
1107 * local_dof_indices_neighbor);
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 *
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 *
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,
1874 * coeffs_re += coeffs_tmp;
1875 *
1876 * coeffs_tmp = 0;
1877 * solver.solve(local_matrix_lift,
1878 * coeffs_tmp,
1879 * local_rhs_be,
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 adjecent 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,
1989 * solver.solve(local_matrix_lift,
1990 * coeffs_be,
1991 * local_rhs_be,
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 > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1064
void reinit(const Triangulation< dim, spacedim > &tria)
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
Definition: fe_values.h:3963
const Quadrature< dim > quadrature
Definition: fe_values.h:4154
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Definition: fe_dgq.h:111
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)
void vmult(Vector< double > &dst, const Vector< double > &src) const
Definition: tensor.h:503
Definition: vector.h:109
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Point< 2 > first
Definition: grid_out.cc:4603
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
Definition: exceptions.h:1473
void write_vtk(std::ostream &out) const
typename ActiveSelector::face_iterator face_iterator
Definition: dof_handler.h:484
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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, SparsityPatternType &sparsity_pattern)
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:190
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)
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:13734
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)