527 *
double return_value = 0.0;
531 * return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
532 * +24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0])) +
533 * 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
534 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]);
538 * return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) *
539 * p[2] * (1.0 - p[2])) +
540 * 24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
541 * p[2] * (1.0 - p[2])) +
542 * 24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
543 * p[1] * (1.0 - p[1])) +
544 * 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
545 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
546 * Utilities::fixed_power<2>(p[2] * (1.0 - p[2])) +
547 * 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
548 * (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
549 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
550 * 2.0 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
551 * (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
552 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
557 *
return return_value;
564 * This
class implement the manufactured (exact) solution @f$u@f$. To compute the
565 * errors, we need the
value of @f$u@f$ as well as its
gradient and its Hessian.
569 *
class ExactSolution :
public Function<dim>
577 *
const unsigned int component = 0)
const override;
581 *
const unsigned int component = 0)
const override;
585 *
const unsigned int component = 0)
const override;
591 *
double ExactSolution<dim>::value(
const Point<dim> &p,
592 *
const unsigned int )
const
594 *
double return_value = 0.0;
599 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
603 * return_value = Utilities::fixed_power<2>(
604 * p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
609 *
return return_value;
616 * ExactSolution<dim>::gradient(
const Point<dim> &p,
617 *
const unsigned int )
const
623 * return_gradient[0] =
624 * (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
625 * 4.0 * Utilities::fixed_power<3>(p[0])) *
626 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
627 * return_gradient[1] =
628 * (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
629 * 4.0 * Utilities::fixed_power<3>(p[1])) *
630 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
634 * return_gradient[0] =
635 * (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
636 * 4.0 * Utilities::fixed_power<3>(p[0])) *
637 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
638 * return_gradient[1] =
639 * (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
640 * 4.0 * Utilities::fixed_power<3>(p[1])) *
641 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
642 * return_gradient[2] =
643 * (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
644 * 4.0 * Utilities::fixed_power<3>(p[2])) *
645 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
650 *
return return_gradient;
657 * ExactSolution<dim>::hessian(
const Point<dim> &p,
658 *
const unsigned int )
const
664 * return_hessian[0][0] = (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
665 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
666 * return_hessian[0][1] =
667 * (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
668 * 4.0 * Utilities::fixed_power<3>(p[0])) *
669 * (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
670 * 4.0 * Utilities::fixed_power<3>(p[1]));
671 * return_hessian[1][1] = (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
672 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
676 * return_hessian[0][0] =
677 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
678 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
679 * return_hessian[0][1] =
680 * (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
681 * 4.0 * Utilities::fixed_power<3>(p[0])) *
682 * (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
683 * 4.0 * Utilities::fixed_power<3>(p[1])) *
684 * Utilities::fixed_power<2>(p[2] * (1.0 - p[2]));
685 * return_hessian[0][2] =
686 * (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
687 * 4.0 * Utilities::fixed_power<3>(p[0])) *
688 * (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
689 * 4.0 * Utilities::fixed_power<3>(p[2])) *
690 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
691 * return_hessian[1][1] =
692 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
693 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
694 * return_hessian[1][2] =
695 * (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
696 * 4.0 * Utilities::fixed_power<3>(p[1])) *
697 * (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
698 * 4.0 * Utilities::fixed_power<3>(p[2])) *
699 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
700 * return_hessian[2][2] =
701 * (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
702 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
707 *
return return_hessian;
715 * <a name=
"step_82-ImplementationofthecodeBiLaplacianLDGLiftcodeclass"></a>
716 * <h3>Implementation of the <code>BiLaplacianLDGLift</code>
class</h3>
721 * <a name=
"step_82-BiLaplacianLDGLiftBiLaplacianLDGLift"></a>
722 * <h4>BiLaplacianLDGLift::BiLaplacianLDGLift</h4>
726 * In the constructor, we set the polynomial degree of the two finite element
727 * spaces, we associate the corresponding DoF handlers to the
triangulation,
728 * and we set the two penalty coefficients.
732 * BiLaplacianLDGLift<dim>::BiLaplacianLDGLift(
const unsigned int n_refinements,
733 *
const unsigned int fe_degree,
734 *
const double penalty_jump_grad,
735 *
const double penalty_jump_val)
736 * : n_refinements(n_refinements)
740 * , penalty_jump_grad(penalty_jump_grad)
741 * , penalty_jump_val(penalty_jump_val)
749 * <a name=
"step_82-BiLaplacianLDGLiftmake_grid"></a>
750 * <h4>BiLaplacianLDGLift::make_grid</h4>
754 * To build a mesh
for @f$\Omega=(0,1)^d@f$, we simply call the function
755 * <code>GridGenerator::hyper_cube</code> and then refine it
using
756 * <code>refine_global</code>. The number of refinements is hard-coded
761 *
void BiLaplacianLDGLift<dim>::make_grid()
763 * std::cout <<
"Building the mesh............." << std::endl;
769 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
778 * <a name=
"step_82-BiLaplacianLDGLiftsetup_system"></a>
779 * <h4>BiLaplacianLDGLift::setup_system</h4>
783 * In the following function, we set up the degrees of freedom, the sparsity
784 * pattern, the size of the
matrix @f$A@f$, and the size of the solution and
785 * right-hand side vectors
786 * @f$\boldsymbol{U}@f$ and @f$\boldsymbol{
F}@f$. For the sparsity pattern, we cannot
788 * (as we would
do for instance
for the SIPG method) because we need to take
789 * into account the interactions of a neighboring cell with another
790 * neighboring cell as described in the introduction. The extended sparsity
791 * pattern is built by iterating over all the active cells. For the current
792 * cell, we collect all its degrees of freedom as well as the degrees of
793 * freedom of all its neighboring cells, and then couple everything with
798 *
void BiLaplacianLDGLift<dim>::setup_system()
800 * dof_handler.distribute_dofs(fe);
802 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
807 *
const auto dofs_per_cell = fe.dofs_per_cell;
809 *
for (
const auto &cell : dof_handler.active_cell_iterators())
812 * cell->get_dof_indices(dofs);
814 *
for (
unsigned int f = 0; f < cell->n_faces(); ++f)
815 *
if (!cell->face(f)->at_boundary())
817 *
const auto neighbor_cell = cell->neighbor(f);
819 * std::vector<types::global_dof_index> tmp(dofs_per_cell);
820 * neighbor_cell->get_dof_indices(tmp);
822 * dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp));
825 *
for (
const auto i : dofs)
826 * for (const auto j : dofs)
833 * sparsity_pattern.copy_from(dsp);
836 *
matrix.reinit(sparsity_pattern);
837 * rhs.reinit(dof_handler.n_dofs());
839 * solution.reinit(dof_handler.n_dofs());
843 * At the
end of the function, we output
this sparsity pattern as
844 * a scalable vector graphic. You can visualize it by loading
this
845 * file in most web browsers:
848 * std::ofstream out(
"sparsity-pattern.svg");
849 * sparsity_pattern.print_svg(out);
857 * <a name=
"step_82-BiLaplacianLDGLiftassemble_system"></a>
858 * <h4>BiLaplacianLDGLift::assemble_system</h4>
862 * This function simply calls the two
functions responsible
863 *
for the assembly of the
matrix and the right-hand side.
867 *
void BiLaplacianLDGLift<dim>::assemble_system()
869 * std::cout <<
"Assembling the system............." << std::endl;
874 * std::cout <<
"Done. " << std::endl;
882 * <a name=
"step_82-BiLaplacianLDGLiftassemble_matrix"></a>
883 * <h4>BiLaplacianLDGLift::assemble_matrix</h4>
887 * This function assembles the
matrix @f$A@f$ whose entries are defined
888 * by @f$A_{ij}=A_h(\varphi_j,\varphi_i)@f$ which involves the product of
889 * discrete Hessians and the penalty terms.
893 *
void BiLaplacianLDGLift<dim>::assemble_matrix()
898 *
const QGauss<dim - 1> quad_face(fe.degree + 1);
900 *
const unsigned int n_q_points = quad.size();
901 *
const unsigned int n_q_points_face = quad_face.size();
911 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
913 * std::vector<types::global_dof_index> local_dof_indices(n_dofs);
914 * std::vector<types::global_dof_index> local_dof_indices_neighbor(n_dofs);
915 * std::vector<types::global_dof_index> local_dof_indices_neighbor_2(n_dofs);
919 * As indicated in the introduction, the following matrices are used
for
920 * the contributions of the products of the discrete Hessians.
938 * The following matrices are used
for the contributions of the two
950 * std::vector<std::vector<Tensor<2, dim>>> discrete_hessians(
952 * std::vector<std::vector<std::vector<Tensor<2, dim>>>>
954 * discrete_hessians);
956 *
for (
const auto &cell : dof_handler.active_cell_iterators())
959 * cell->get_dof_indices(local_dof_indices);
963 * We now compute all the discrete Hessians that are not vanishing
964 * on the current cell, i.e., the discrete Hessian of all the basis
965 *
functions with support on the current cell or on one of its
969 * compute_discrete_hessians(cell,
971 * discrete_hessians_neigh);
975 * First, we compute and add the interactions of the degrees of freedom
976 * of the current cell.
979 * stiffness_matrix_cc = 0;
980 *
for (
unsigned int q = 0; q < n_q_points; ++q)
982 *
const double dx = fe_values.JxW(q);
984 *
for (
unsigned int i = 0; i < n_dofs; ++i)
985 *
for (
unsigned int j = 0; j < n_dofs; ++j)
990 * stiffness_matrix_cc(i, j) += scalar_product(H_j, H_i) *
dx;
994 *
for (
unsigned int i = 0; i < n_dofs; ++i)
995 *
for (
unsigned int j = 0; j < n_dofs; ++j)
997 *
matrix(local_dof_indices[i], local_dof_indices[j]) +=
998 * stiffness_matrix_cc(i, j);
1003 * Next, we compute and add the interactions of the degrees of freedom
1004 * of the current cell with those of its neighbors. Note that the
1005 * interactions of the degrees of freedom of a neighbor with those of
1006 * the same neighbor are included here.
1009 *
for (
unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1012 * cell->face(face_no);
1014 *
const bool at_boundary = face->at_boundary();
1019 * There is
nothing to be done
if boundary face (the liftings of
1020 * the Dirichlet BCs are accounted
for in the assembly of the
1021 * RHS; in fact,
nothing to be done in
this program since we
1022 * prescribe homogeneous BCs).
1029 * neighbor_cell = cell->neighbor(face_no);
1030 * neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
1032 * stiffness_matrix_cn = 0;
1033 * stiffness_matrix_nc = 0;
1034 * stiffness_matrix_nn = 0;
1035 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1037 *
const double dx = fe_values.JxW(q);
1039 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1041 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1047 * discrete_hessians_neigh[face_no][i][q];
1049 * discrete_hessians_neigh[face_no][j][q];
1051 * stiffness_matrix_cn(i, j) +=
1052 * scalar_product(H_j_neigh, H_i) *
dx;
1053 * stiffness_matrix_nc(i, j) +=
1054 * scalar_product(H_j, H_i_neigh) *
dx;
1055 * stiffness_matrix_nn(i, j) +=
1056 * scalar_product(H_j_neigh, H_i_neigh) *
dx;
1061 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1063 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1065 *
matrix(local_dof_indices[i],
1066 * local_dof_indices_neighbor[j]) +=
1067 * stiffness_matrix_cn(i, j);
1068 *
matrix(local_dof_indices_neighbor[i],
1069 * local_dof_indices[j]) +=
1070 * stiffness_matrix_nc(i, j);
1071 *
matrix(local_dof_indices_neighbor[i],
1072 * local_dof_indices_neighbor[j]) +=
1073 * stiffness_matrix_nn(i, j);
1082 * We now compute and add the interactions of the degrees of freedom of
1083 * a neighboring cells with those of another neighboring cell (
this is
1084 * where we need the extended sparsity pattern).
1087 *
for (
unsigned int face_no = 0; face_no < cell->n_faces() - 1; ++face_no)
1090 * cell->face(face_no);
1092 *
const bool at_boundary = face->at_boundary();
1097 * Dirichlet BCs are accounted
for in the assembly of the RHS;
1098 * in fact,
nothing to be done in
this program since we
1099 * prescribe homogeneous BCs)
1108 *
for (
unsigned int face_no_2 = face_no + 1;
1109 * face_no_2 < cell->n_faces();
1113 * cell->face(face_no_2);
1115 *
const bool at_boundary_2 = face_2->at_boundary();
1116 *
if (!at_boundary_2)
1119 * neighbor_cell = cell->neighbor(face_no);
1120 * neighbor_cell->get_dof_indices(
1121 * local_dof_indices_neighbor);
1123 * neighbor_cell_2 = cell->neighbor(face_no_2);
1124 * neighbor_cell_2->get_dof_indices(
1125 * local_dof_indices_neighbor_2);
1127 * stiffness_matrix_n1n2 = 0;
1128 * stiffness_matrix_n2n1 = 0;
1130 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1132 *
const double dx = fe_values.JxW(q);
1134 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1135 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1138 * discrete_hessians_neigh[face_no][i][q];
1140 * discrete_hessians_neigh[face_no][j][q];
1143 * discrete_hessians_neigh[face_no_2][i][q];
1145 * discrete_hessians_neigh[face_no_2][j][q];
1147 * stiffness_matrix_n1n2(i, j) +=
1148 * scalar_product(H_j_neigh2, H_i_neigh) *
dx;
1149 * stiffness_matrix_n2n1(i, j) +=
1150 * scalar_product(H_j_neigh, H_i_neigh2) *
dx;
1154 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1155 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1157 *
matrix(local_dof_indices_neighbor[i],
1158 * local_dof_indices_neighbor_2[j]) +=
1159 * stiffness_matrix_n1n2(i, j);
1160 *
matrix(local_dof_indices_neighbor_2[i],
1161 * local_dof_indices_neighbor[j]) +=
1162 * stiffness_matrix_n2n1(i, j);
1172 * Finally, we compute and add the two penalty terms.
1175 *
for (
unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1178 * cell->face(face_no);
1180 *
const double mesh_inv = 1.0 / face->diameter();
1181 *
const double mesh3_inv =
1182 * 1.0 / Utilities::fixed_power<3>(face->diameter());
1184 * fe_face.reinit(cell, face_no);
1188 *
const bool at_boundary = face->at_boundary();
1191 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1193 *
const double dx = fe_face.JxW(q);
1195 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1196 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1198 * ip_matrix_cc(i, j) += penalty_jump_grad * mesh_inv *
1199 * fe_face.shape_grad(j, q) *
1200 * fe_face.shape_grad(i, q) *
dx;
1201 * ip_matrix_cc(i, j) += penalty_jump_val * mesh3_inv *
1202 * fe_face.shape_value(j, q) *
1203 * fe_face.shape_value(i, q) *
dx;
1211 * neighbor_cell = cell->neighbor(face_no);
1212 *
const unsigned int face_no_neighbor =
1213 * cell->neighbor_of_neighbor(face_no);
1217 * In the next step, we need to have a global way to compare the
1218 * cells in order to not calculate the same jump term twice:
1221 *
if (neighbor_cell->id() < cell->id())
1225 * fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1226 * neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
1232 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1234 *
const double dx = fe_face.JxW(q);
1236 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1238 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1240 * ip_matrix_cc(i, j) +=
1241 * penalty_jump_grad * mesh_inv *
1242 * fe_face.shape_grad(j, q) *
1243 * fe_face.shape_grad(i, q) *
dx;
1244 * ip_matrix_cc(i, j) +=
1245 * penalty_jump_val * mesh3_inv *
1246 * fe_face.shape_value(j, q) *
1247 * fe_face.shape_value(i, q) *
dx;
1249 * ip_matrix_cn(i, j) -=
1250 * penalty_jump_grad * mesh_inv *
1251 * fe_face_neighbor.shape_grad(j, q) *
1252 * fe_face.shape_grad(i, q) *
dx;
1253 * ip_matrix_cn(i, j) -=
1254 * penalty_jump_val * mesh3_inv *
1255 * fe_face_neighbor.shape_value(j, q) *
1256 * fe_face.shape_value(i, q) *
dx;
1258 * ip_matrix_nc(i, j) -=
1259 * penalty_jump_grad * mesh_inv *
1260 * fe_face.shape_grad(j, q) *
1261 * fe_face_neighbor.shape_grad(i, q) *
dx;
1262 * ip_matrix_nc(i, j) -=
1263 * penalty_jump_val * mesh3_inv *
1264 * fe_face.shape_value(j, q) *
1265 * fe_face_neighbor.shape_value(i, q) *
dx;
1267 * ip_matrix_nn(i, j) +=
1268 * penalty_jump_grad * mesh_inv *
1269 * fe_face_neighbor.shape_grad(j, q) *
1270 * fe_face_neighbor.shape_grad(i, q) *
dx;
1271 * ip_matrix_nn(i, j) +=
1272 * penalty_jump_val * mesh3_inv *
1273 * fe_face_neighbor.shape_value(j, q) *
1274 * fe_face_neighbor.shape_value(i, q) *
dx;
1282 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1284 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1286 *
matrix(local_dof_indices[i], local_dof_indices[j]) +=
1287 * ip_matrix_cc(i, j);
1293 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1295 *
for (
unsigned int j = 0; j < n_dofs; ++j)
1297 *
matrix(local_dof_indices[i],
1298 * local_dof_indices_neighbor[j]) +=
1299 * ip_matrix_cn(i, j);
1300 *
matrix(local_dof_indices_neighbor[i],
1301 * local_dof_indices[j]) += ip_matrix_nc(i, j);
1302 *
matrix(local_dof_indices_neighbor[i],
1303 * local_dof_indices_neighbor[j]) +=
1304 * ip_matrix_nn(i, j);
1318 * <a name=
"step_82-BiLaplacianLDGLiftassemble_rhs"></a>
1319 * <h4>BiLaplacianLDGLift::assemble_rhs</h4>
1323 * This function
assemble the right-hand side of the system. Since we consider
1324 * homogeneous Dirichlet boundary conditions, imposed weakly in the bilinear
1325 * form
using the Nitsche approach, it only involves the contribution of the
1326 * forcing term @f$\int_{\Omega}fv_h@f$.
1329 *
template <
int dim>
1330 *
void BiLaplacianLDGLift<dim>::assemble_rhs()
1338 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
1339 *
const unsigned int n_quad_pts = quad.size();
1341 *
const RightHandSide<dim> right_hand_side;
1344 * std::vector<types::global_dof_index> local_dof_indices(n_dofs);
1346 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1348 * fe_values.
reinit(cell);
1349 * cell->get_dof_indices(local_dof_indices);
1352 *
for (
unsigned int q = 0; q < n_quad_pts; ++q)
1354 *
const double dx = fe_values.JxW(q);
1356 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1359 * right_hand_side.value(fe_values.quadrature_point(q)) *
1360 * fe_values.shape_value(i, q) *
dx;
1364 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1365 * rhs(local_dof_indices[i]) += local_rhs(i);
1374 * <a name=
"step_82-BiLaplacianLDGLiftsolve"></a>
1375 * <h4>BiLaplacianLDGLift::solve</h4>
1379 * To solve the linear system @f$A\boldsymbol{U}=\boldsymbol{
F}@f$,
1380 * we proceed as in @ref step_74
"step-74" and use a direct method. We could
1381 * as well use an iterative method,
for instance the conjugate
1382 *
gradient method as in @ref step_3
"step-3".
1385 *
template <
int dim>
1386 *
void BiLaplacianLDGLift<dim>::solve()
1390 * A_direct.vmult(solution, rhs);
1398 * <a name=
"step_82-BiLaplacianLDGLiftcompute_errors"></a>
1399 * <h4>BiLaplacianLDGLift::compute_errors</h4>
1403 * This function computes the discrete @f$H^2@f$, @f$H^1@f$ and @f$L^2@f$ norms of
1404 * the error @f$u-u_h@f$, where @f$u@f$ is the exact solution and @f$u_h@f$ is
1405 * the
approximate solution. See the introduction
for the definition
1409 *
template <
int dim>
1410 *
void BiLaplacianLDGLift<dim>::compute_errors()
1412 *
double error_H2 = 0;
1413 *
double error_H1 = 0;
1414 *
double error_L2 = 0;
1417 *
const QGauss<dim - 1> quad_face(fe.degree + 1);
1433 *
const unsigned int n_q_points = quad.size();
1434 *
const unsigned int n_q_points_face = quad_face.size();
1438 * We introduce some variables
for the exact solution...
1441 *
const ExactSolution<dim> u_exact;
1448 * std::vector<double> solution_values_cell(n_q_points);
1449 * std::vector<Tensor<1, dim>> solution_gradients_cell(n_q_points);
1450 * std::vector<Tensor<2, dim>> solution_hessians_cell(n_q_points);
1452 * std::vector<double> solution_values(n_q_points_face);
1453 * std::vector<double> solution_values_neigh(n_q_points_face);
1454 * std::vector<Tensor<1, dim>> solution_gradients(n_q_points_face);
1455 * std::vector<Tensor<1, dim>> solution_gradients_neigh(n_q_points_face);
1457 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1459 * fe_values.
reinit(cell);
1461 * fe_values.get_function_values(solution, solution_values_cell);
1462 * fe_values.get_function_gradients(solution, solution_gradients_cell);
1463 * fe_values.get_function_hessians(solution, solution_hessians_cell);
1467 * We
first add the <i>bulk</i> terms.
1470 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1472 *
const double dx = fe_values.JxW(q);
1474 * error_H2 += (u_exact.hessian(fe_values.quadrature_point(q)) -
1475 * solution_hessians_cell[q])
1478 * error_H1 += (u_exact.gradient(fe_values.quadrature_point(q)) -
1479 * solution_gradients_cell[q])
1482 * error_L2 += Utilities::fixed_power<2>(
1483 * u_exact.value(fe_values.quadrature_point(q)) -
1484 * solution_values_cell[q]) *
1490 * We then add the face contributions.
1493 *
for (
unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1496 * cell->face(face_no);
1498 *
const double mesh_inv = 1.0 / face->diameter();
1499 *
const double mesh3_inv =
1500 * 1.0 / Utilities::fixed_power<3>(face->diameter());
1502 * fe_face.reinit(cell, face_no);
1504 * fe_face.get_function_values(solution, solution_values);
1505 * fe_face.get_function_gradients(solution, solution_gradients);
1507 *
const bool at_boundary = face->at_boundary();
1510 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1512 *
const double dx = fe_face.JxW(q);
1513 *
const double u_exact_q =
1514 * u_exact.value(fe_face.quadrature_point(q));
1516 * u_exact.gradient(fe_face.quadrature_point(q));
1520 * (u_exact_grad_q - solution_gradients[q]).
norm_square() *
1522 * error_H2 += mesh3_inv *
1523 * Utilities::fixed_power<2>(u_exact_q -
1524 * solution_values[q]) *
1526 * error_H1 += mesh_inv *
1527 * Utilities::fixed_power<2>(u_exact_q -
1528 * solution_values[q]) *
1536 * neighbor_cell = cell->neighbor(face_no);
1537 *
const unsigned int face_no_neighbor =
1538 * cell->neighbor_of_neighbor(face_no);
1542 * In the next step, we need to have a global way to compare the
1543 * cells in order to not calculate the same jump term twice:
1546 *
if (neighbor_cell->id() < cell->id())
1550 * fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1552 * fe_face.get_function_values(solution, solution_values);
1553 * fe_face_neighbor.get_function_values(solution,
1554 * solution_values_neigh);
1555 * fe_face.get_function_gradients(solution,
1556 * solution_gradients);
1557 * fe_face_neighbor.get_function_gradients(
1558 * solution, solution_gradients_neigh);
1560 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1562 *
const double dx = fe_face.JxW(q);
1566 * To compute the jump term, we use the fact that
1567 * @f$\jump{u}=0@f$ and
1568 * @f$\jump{\nabla u}=\mathbf{0}@f$ since @f$u\in
1574 * (solution_gradients_neigh[q] - solution_gradients[q])
1579 * Utilities::fixed_power<2>(solution_values_neigh[q] -
1580 * solution_values[q]) *
1584 * Utilities::fixed_power<2>(solution_values_neigh[q] -
1585 * solution_values[q]) *
1600 * std::cout <<
"DG H2 norm of the error: " << error_H2 << std::endl;
1601 * std::cout <<
"DG H1 norm of the error: " << error_H1 << std::endl;
1602 * std::cout <<
" L2 norm of the error: " << error_L2 << std::endl;
1610 * <a name=
"step_82-BiLaplacianLDGLiftoutput_results"></a>
1611 * <h4>BiLaplacianLDGLift::output_results</h4>
1615 * This function, which writes the solution to a
vtk file,
1616 * is copied from @ref step_3
"step-3".
1619 *
template <
int dim>
1620 *
void BiLaplacianLDGLift<dim>::output_results() const
1624 * data_out.add_data_vector(solution,
"solution");
1625 * data_out.build_patches();
1627 * std::ofstream output(
"solution.vtk");
1628 * data_out.write_vtk(output);
1636 * <a name=
"step_82-BiLaplacianLDGLiftassemble_local_matrix"></a>
1637 * <h4>BiLaplacianLDGLift::assemble_local_matrix</h4>
1641 * As already mentioned above,
this function is used to
assemble
1642 * the (local) mass matrices needed
for the computations of the
1643 * lifting terms. We reiterate that only the basis
functions with
1644 * support on the current cell are considered.
1647 *
template <
int dim>
1648 *
void BiLaplacianLDGLift<dim>::assemble_local_matrix(
1650 *
const unsigned int n_q_points,
1655 *
const unsigned int n_dofs = fe_values_lift.dofs_per_cell;
1658 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1660 *
const double dx = fe_values_lift.JxW(q);
1662 *
for (
unsigned int m = 0; m < n_dofs; ++m)
1663 *
for (
unsigned int n = 0; n < n_dofs; ++n)
1665 * local_matrix(m, n) +=
1666 * scalar_product(fe_values_lift[tau_ext].
value(n, q),
1667 * fe_values_lift[tau_ext].
value(m, q)) *
1678 * <a name=
"step_82-BiLaplacianLDGLiftcompute_discrete_hessians"></a>
1679 * <h4>BiLaplacianLDGLift::compute_discrete_hessians</h4>
1683 * This function is the main novelty of
this program. It computes
1684 * the discrete Hessian @f$H_h(\varphi)@f$
for all the basis
functions
1685 * @f$\varphi@f$ of @f$\mathbb{V}_h@f$ supported on the current cell and
1686 * those supported on a neighboring cell. The
first argument
1687 * indicates the current cell (referring to the global
DoFHandler
1688 *
object),
while the other two arguments are output variables that
1689 * are filled by
this function.
1693 * In the following, we need to evaluate finite element shape
1694 *
functions for the `fe_lift` finite element on the current
1695 * cell. Like
for example in @ref step_61
"step-61",
this "lift" space is defined
1696 * on every cell individually; as a consequence, there is no global
1697 *
DoFHandler associated with
this because we simply have no need
1698 *
for such a
DoFHandler. That leaves the question of what we should
1700 * them to evaluate shape
functions of `fe_lift` on a concrete
1701 * cell. If we simply provide the
first argument to
this function,
1703 * that the given `cell` belongs to a
DoFHandler that has a
1704 * different finite element associated with it than the `fe_lift`
1705 *
object we want to evaluate. Fortunately, there is a relatively
1707 * points into a
triangulation -- the same cell, but not associated
1708 * with a
DoFHandler, and consequently no finite element space. In
1709 * that case,
FEValues::reinit() will skip the check that would
1710 * otherwise lead to an error message. All we have to do is to convert
1712 * see the
first couple of lines of the function below to see how
1716 * template <
int dim>
1717 *
void BiLaplacianLDGLift<dim>::compute_discrete_hessians(
1718 * const typename
DoFHandler<dim>::active_cell_iterator &cell,
1719 *
std::vector<
std::vector<
Tensor<2, dim>>> &discrete_hessians,
1721 * &discrete_hessians_neigh)
1727 *
const QGauss<dim - 1> quad_face(
fe.degree + 1);
1729 *
const unsigned int n_q_points = quad.size();
1730 *
const unsigned int n_q_points_face = quad_face.size();
1734 * The information we need from the basis
functions of
1735 * @f$\mathbb{V}_h@f$: <code>fe_values</code> is needed to add
1736 * the broken Hessian part of the discrete Hessian,
while
1737 * <code>fe_face</code> and <code>fe_face_neighbor</code>
1738 * are used to compute the right-hand sides
for the local
1750 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
1754 * The information needed from the basis
functions
1755 * of the finite element space
for the lifting terms:
1756 * <code>fe_values_lift</code> is used
for the (local)
1757 * mass
matrix (see @f$\boldsymbol{M}_c@f$ in the introduction),
1758 * while <code>fe_face_lift</code> is used to compute the
1759 * right-hand sides (see @f$\boldsymbol{G}_c@f$
for @f$b_e@f$).
1771 *
const unsigned int n_dofs_lift = fe_values_lift.dofs_per_cell;
1774 *
Vector<double> local_rhs_re(n_dofs_lift), local_rhs_be(n_dofs_lift),
1775 * coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift), coeffs_tmp(n_dofs_lift);
1780 *
double factor_avg;
1782 * fe_values.reinit(cell);
1783 * fe_values_lift.reinit(cell_lift);
1787 * We start by assembling the (local) mass
matrix used
for the computation
1788 * of the lifting terms @f$r_e@f$ and @f$b_e@f$.
1791 * assemble_local_matrix(fe_values_lift, n_q_points, local_matrix_lift);
1793 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1794 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1796 * discrete_hessians[i][q] = 0;
1798 *
for (
unsigned int face_no = 0;
1799 * face_no < discrete_hessians_neigh.size();
1802 * discrete_hessians_neigh[face_no][i][q] = 0;
1809 * @f$x_q@f$ of <code>cell</code>
for each basis function supported on
1810 * <code>cell</code>, namely we fill-in the variable
1811 * <code>discrete_hessians[i][q]</code>. For the lifting terms, we need to
1812 * add the contribution of all the faces of <code>cell</code>.
1815 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1820 *
for (
unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1823 * cell->face(face_no);
1825 *
const bool at_boundary = face->at_boundary();
1829 * Recall that by convention, the average of a function across a
1830 * boundary face @f$e@f$ reduces to the
trace of the function on the
1831 * only element adjacent to @f$e@f$, namely there is no factor
1832 * @f$\frac{1}{2}@f$. We distinguish between the two cases (the current
1833 * face lies in the interior or on the boundary of the domain)
using
1834 * the variable <code>factor_avg</code>.
1843 * fe_face.reinit(cell, face_no);
1844 * fe_face_lift.reinit(cell_lift, face_no);
1847 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1849 *
const double dx = fe_face_lift.JxW(q);
1853 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1855 * local_rhs_re(m) +=
1857 * (fe_face_lift[tau_ext].
value(m, q) * normal) *
1858 * fe_face.shape_grad(i, q) *
dx;
1864 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
1865 * introduced in the comments about the implementation of the
1866 * lifting @f$b_e@f$ in the
case
1867 * @f$\varphi=\varphi^c@f$.
1871 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1873 *
const double dx = fe_face_lift.JxW(q);
1877 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1879 * local_rhs_be(m) += factor_avg *
1880 * fe_face_lift[tau_ext].divergence(m, q) *
1881 * normal * fe_face.shape_value(i, q) *
dx;
1886 * solver.solve(local_matrix_lift,
1890 * coeffs_re += coeffs_tmp;
1893 * solver.solve(local_matrix_lift,
1897 * coeffs_be += coeffs_tmp;
1901 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1903 * discrete_hessians[i][q] += fe_values.shape_hessian(i, q);
1905 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1907 * discrete_hessians[i][q] -=
1908 * coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
1911 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1913 * discrete_hessians[i][q] +=
1914 * coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
1924 * @f$x_q@f$ of <code>cell</code>
for each basis function supported on a
1925 * neighboring <code>neighbor_cell</code> of <code>cell</code>, namely we
1926 * fill-in the variable <code>discrete_hessians_neigh[face_no][i][q]</code>.
1927 * For the lifting terms, we only need to add the contribution of the
1928 * face adjacent to <code>cell</code> and <code>neighbor_cell</code>.
1931 *
for (
unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
1934 * cell->face(face_no);
1936 *
const bool at_boundary = face->at_boundary();
1942 * For non-homogeneous Dirichlet BCs, we would need to
1943 * compute the lifting of the prescribed BC (see the
1944 *
"Possible Extensions" section
for more details).
1951 * neighbor_cell = cell->neighbor(face_no);
1952 *
const unsigned int face_no_neighbor =
1953 * cell->neighbor_of_neighbor(face_no);
1954 * fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
1956 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1961 * fe_face_lift.reinit(cell_lift, face_no);
1964 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1966 *
const double dx = fe_face_lift.JxW(q);
1968 * fe_face_neighbor.normal_vector(q);
1970 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1972 * local_rhs_re(m) +=
1973 * 0.5 * (fe_face_lift[tau_ext].value(m, q) * normal) *
1974 * fe_face_neighbor.shape_grad(i, q) *
dx;
1980 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
1981 * introduced in the comments about the implementation of the
1982 * lifting @f$b_e@f$ in the
case
1983 * @f$\varphi=\varphi^n@f$.
1987 *
for (
unsigned int q = 0; q < n_q_points_face; ++q)
1989 *
const double dx = fe_face_lift.JxW(q);
1991 * fe_face_neighbor.normal_vector(q);
1993 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
1995 * local_rhs_be(m) +=
1996 * 0.5 * fe_face_lift[tau_ext].divergence(m, q) *
1997 * normal * fe_face_neighbor.shape_value(i, q) *
dx;
2001 * solver.solve(local_matrix_lift,
2005 * solver.solve(local_matrix_lift,
2010 *
for (
unsigned int q = 0; q < n_q_points; ++q)
2012 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
2014 * discrete_hessians_neigh[face_no][i][q] -=
2015 * coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
2018 *
for (
unsigned int m = 0; m < n_dofs_lift; ++m)
2020 * discrete_hessians_neigh[face_no][i][q] +=
2021 * coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
2035 * <a name=
"step_82-BiLaplacianLDGLiftrun"></a>
2036 * <h4>BiLaplacianLDGLift::run</h4>
2039 *
template <
int dim>
2040 *
void BiLaplacianLDGLift<dim>::run()
2045 * assemble_system();
2060 * <a name=
"step_82-Thecodemaincodefunction"></a>
2061 * <h3>The <code>main</code> function</h3>
2065 * This is the <code>main</code> function. We define here the number of mesh
2066 * refinements, the polynomial degree
for the two finite element spaces
2067 * (
for the solution and the two liftings) and the two penalty coefficients.
2068 * We can also change the dimension to run the code in 3
d.
2075 *
const unsigned int n_ref = 3;
2077 *
const unsigned int degree =
2080 *
const double penalty_grad =
2082 *
const double penalty_val =
2085 * Step82::BiLaplacianLDGLift<2> problem(n_ref,
2092 *
catch (std::exception &exc)
2094 * std::cerr << std::endl
2096 * <<
"----------------------------------------------------"
2098 * std::cerr <<
"Exception on processing: " << std::endl
2099 * << exc.what() << std::endl
2100 * <<
"Aborting!" << std::endl
2101 * <<
"----------------------------------------------------"
2107 * std::cerr << std::endl
2109 * <<
"----------------------------------------------------"
2111 * std::cerr <<
"Unknown exception!" << std::endl
2112 * <<
"Aborting!" << std::endl
2113 * <<
"----------------------------------------------------"
2121<a name=
"step_82-Results"></a><h1>Results</h1>
2125When running the program, the sparsity pattern is written to an
svg file, the solution is written to a
vtk file, and some results are printed to the console. With the current setup, the output should read
2129Number of active cells: 64
2130Number of degrees of freedom: 576
2131Assembling the system.............
2133DG H2
norm of the error: 0.0151063
2134DG H1
norm of the error: 0.000399747
2135 L2 norm of the error: 5.33856e-05
2139This corresponds to the bi-Laplacian problem with the manufactured solution mentioned above
for @f$d=2@f$, 3 refinements of the mesh, degree @f$k=2@f$, and @f$\gamma_0=\gamma_1=1@f$
for the penalty coefficients. By changing the number of refinements, we get the following results:
2141<table align=
"center" class=
"doxtable">
2154 <td align=
"center">1</td>
2155 <td align=
"right">4</td>
2156 <td align=
"right">36</td>
2157 <td align=
"center">5.651e-02</td>
2158 <td align=
"center">--</td>
2159 <td align=
"center">3.366e-03</td>
2160 <td align=
"center">--</td>
2161 <td align=
"center">3.473e-04</td>
2162 <td align=
"center">--</td>
2165 <td align=
"center">2</td>
2166 <td align=
"right">16</td>
2167 <td align=
"right">144</td>
2168 <td align=
"center">3.095e-02</td>
2169 <td align=
"center">0.87</td>
2170 <td align=
"center">1.284e-03</td>
2171 <td align=
"center">1.39</td>
2172 <td align=
"center">1.369e-04</td>
2173 <td align=
"center">1.34</td>
2176 <td align=
"center">3</td>
2177 <td align=
"right">64</td>
2178 <td align=
"right">576</td>
2179 <td align=
"center">1.511e-02</td>
2180 <td align=
"center">1.03</td>
2181 <td align=
"center">3.997e-04</td>
2182 <td align=
"center">1.68</td>
2183 <td align=
"center">5.339e-05</td>
2184 <td align=
"center">1.36</td>
2187 <td align=
"center">4</td>
2188 <td align=
"right">256</td>
2189 <td align=
"right">2304</td>
2190 <td align=
"center">7.353e-03</td>
2191 <td align=
"center">1.04</td>
2192 <td align=
"center">1.129e-04</td>
2193 <td align=
"center">1.82</td>
2194 <td align=
"center">1.691e-05</td>
2195 <td align=
"center">1.66</td>
2198 <td align=
"center">5</td>
2199 <td align=
"right">1024</td>
2200 <td align=
"right">9216</td>
2201 <td align=
"center">3.609e-03</td>
2202 <td align=
"center">1.03</td>
2203 <td align=
"center">3.024e-05</td>
2204 <td align=
"center">1.90</td>
2205 <td align=
"center">4.789e-06</td>
2206 <td align=
"center">1.82</td>
2209 <td align=
"center">6</td>
2210 <td align=
"right">4096</td>
2211 <td align=
"right">36864</td>
2212 <td align=
"center">1.785e-03</td>
2213 <td align=
"center">1.02</td>
2214 <td align=
"center">7.850e-06</td>
2215 <td align=
"center">1.95</td>
2216 <td align=
"center">1.277e-06</td>
2217 <td align=
"center">1.91</td>
2221This matches the expected optimal convergence rates
for the @f$H^2@f$ and
2222@f$H^1@f$ norms, but is sub-optimal
for the @f$L_2@f$
norm. Incidentally,
this
2223also matches the results seen in @ref step_47
"step-47" when
using polynomial degree
2226Indeed, just like in @ref step_47
"step-47", we can regain the optimal convergence
2227order
if we set the polynomial degree of the finite elements to @f$k=3@f$
2228or higher. Here are the
numbers for @f$k=3@f$:
2230<table align=
"center" class=
"doxtable">
2231 <tr> <th> n_ref </th> <th>
n_cells </th> <th> n_dofs </th> <th> error H2 </th> <th> rate </th> <th> error H1 </th> <th> rate </th> <th> error
L2 </th> <th> rate</th> </tr>
2232 <tr> <td> 1 </td> <td> 4 </td> <td> 36 </td> <td> 1.451e-02 </td> <td> -- </td> <td> 5.494e-04 </td> <td> -- </td> <td> 3.035e-05 </td> <td> --</td> </tr>
2233 <tr> <td> 2 </td> <td> 16 </td> <td> 144 </td> <td> 3.565e-03 </td> <td> 2.02 </td> <td> 6.870e-05 </td> <td> 3.00 </td> <td> 2.091e-06 </td> <td> 3.86</td> </tr>
2234 <tr> <td> 3 </td> <td> 64 </td> <td> 576 </td> <td> 8.891e-04 </td> <td> 2.00 </td> <td> 8.584e-06 </td> <td> 3.00 </td> <td> 1.352e-07 </td> <td> 3.95</td> </tr>
2235 <tr> <td> 4 </td> <td> 256 </td> <td> 2304 </td> <td> 2.223e-04 </td> <td> 2.00 </td> <td> 1.073e-06 </td> <td> 3.00 </td> <td> 8.594e-09 </td> <td> 3.98</td> </tr>
2236 <tr> <td> 5 </td> <td> 1024 </td> <td> 9216 </td> <td> 5.560e-05 </td> <td> 2.00 </td> <td> 1.341e-07 </td> <td> 3.00 </td> <td> 5.418e-10 </td> <td> 3.99</td> </tr>
2237 <tr> <td> 6 </td> <td> 4096 </td> <td> 36864 </td> <td> 1.390e-05 </td> <td> 2.00 </td> <td> 1.676e-08 </td> <td> 3.00 </td> <td> 3.245e-11 </td> <td> 4.06</td> </tr>
2241<a name=
"step_82-Possibleextensions"></a><h3>Possible extensions</h3>
2244The code can be easily adapted to deal with the following cases:
2247 <li>Non-homogeneous Dirichlet boundary conditions on (part of) the boundary @f$\partial \Omega@f$ of @f$\Omega@f$.</li>
2248 <li>Hanging-nodes (proceed as in @ref step_14
"step-14" to not visit a sub-face twice when computing the lifting terms in <code>compute_discrete_hessians</code> and the penalty terms in <code>assemble_matrix</code>).</li>
2249 <li>LDG method
for the Poisson problem (use the discrete gradient consisting of the broken gradient and the lifting of the jump of @f$u_h@f$).</li>
2252We give below additional details
for the
first of these points.
2255<a name=
"step_82-NonhomogeneousDirichletboundaryconditions"></a><h4>Non-homogeneous Dirichlet boundary conditions</h4>
2257If we prescribe non-homogeneous Dirichlet conditions, say
2259\nabla u=\mathbf{g} \quad \mbox{and} \quad u=g \qquad \mbox{on } \partial \Omega,
2261then the right-hand side @f$\boldsymbol{
F}@f$ of the linear system needs to be modified as follows
2263F_i:=\int_{\Omega}f\varphi_i-\sum_{e\in\mathcal{
E}_h^
b}\int_{\Omega}r_e(\mathbf{g}):H_h(\varphi_i)+\sum_{
e\in\mathcal{
E}_h^
b}\int_{\Omega}b_e(g):H_h(\varphi_i)+\gamma_1\sum_{
e\in\mathcal{
E}_h^
b}h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0\sum_{e\in\mathcal{
E}_h^
b}h_e^{-3}\int_e g\varphi_i, \qquad 1\leq i \leq N_h.
2265Note that
for any given
index @f$i@f$, many of the terms are zero. Indeed,
for @f$e\in \mathcal{
E}_h^
b@f$ we have @f${\rm supp}\,(r_e(\mathbf{g}))={\rm supp}\,(b_e(g))=K@f$, where @f$K@f$ is the element
for which @f$e\subset\partial K@f$. Therefore, the liftings @f$r_e(\mathbf{g})@f$ and @f$b_e(g)@f$ contribute to @f$F_i@f$ only
if @f$\varphi_i@f$ has support on @f$K@f$ or a neighbor of @f$K@f$. In other words, when integrating on a cell @f$K@f$, we need to add
2267\int_{
K}f\varphi_i+\sum_{e\in\mathcal{
E}_h^
b, e\subset\partial
K}\left[-\int_{
K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{
K}b_e(g):H_h(\varphi_i)+\gamma_1h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0h_e^{-3}\int_e g\varphi_i\right]
2269to @f$F_i@f$
for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on @f$K@f$ and
2271\sum_{e\in\mathcal{
E}_h^
b, e\subset\partial
K}\left[-\int_{
K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{
K}b_e(g):H_h(\varphi_i)\right]
2273to @f$F_i@f$ for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on a neighbor of @f$K@f$.
2276Note that we can easily consider the case where Dirichlet boundary conditions are imposed only on a subset @f$\emptyset\neq\Gamma_D\subset\partial \Omega@f$. In this case, we simply need to replace @f$\mathcal{
E}_h^
b@f$ by @f$\mathcal{
E}_h^D\subset\mathcal{
E}_h^
b@f$ consisting of the faces belonging to @f$\Gamma_D@f$. This also affects the
matrix @f$A@f$ (simply replace @f$\mathcal{
E}_h=\mathcal{
E}_h^0\cup\mathcal{
E}_h^
b@f$ by @f$\mathcal{
E}_h=\mathcal{
E}_h^0\cup\mathcal{
E}_h^D@f$).
2279<a name=
"step_82-PlainProg"></a>
2280<h1> The plain program</h1>
2281@include
"step-82.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
const Quadrature< dim > quadrature
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
void initialize(const SparsityPattern &sparsity_pattern)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)