From 2ca9d4d83dda73ff6c058bebf4cc77ace4026708 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 2 Sep 2021 22:20:24 -0600 Subject: [PATCH] Use the usual style in main(). --- examples/step-82/step-82.cc | 1805 +++++++++++++++++++++-------------- 1 file changed, 1066 insertions(+), 739 deletions(-) diff --git a/examples/step-82/step-82.cc b/examples/step-82/step-82.cc index 6e84beb76e..5071d9967a 100644 --- a/examples/step-82/step-82.cc +++ b/examples/step-82/step-82.cc @@ -1,4 +1,3 @@ - /* --------------------------------------------------------------------- * * Copyright (C) 2021 by the deal.II authors @@ -58,7 +57,6 @@ namespace Step82 { - using namespace dealii; // @sect3{The BiLaplacianLDGLift class template} @@ -71,13 +69,13 @@ namespace Step82 class BiLaplacianLDGLift { public: - - BiLaplacianLDGLift(const unsigned int fe_degree, double penalty_jump_grad, double penalty_jump_val); + BiLaplacianLDGLift(const unsigned int fe_degree, + double penalty_jump_grad, + double penalty_jump_val); void run(); private: - void make_grid(); void setup_system(); void assemble_system(); @@ -85,7 +83,7 @@ namespace Step82 void assemble_rhs(); void solve(); - + void compute_errors(); void output_results() const; @@ -93,24 +91,30 @@ namespace Step82 // is used for the assembly of the (local) mass matrix used to compute the // two lifting terms (see the matrix $\boldsymbol{M}_c$ introduced in // the introduction when describing the computation of $b_e$). The function - // compute_discrete_hessians computes the required discrete Hessians: - // the discrete Hessians of the basis functions with support on the current - // cell (stored in the output variable discrete_hessians) - // and the basis functions with support on a neighbor of the current cell - // (stored in the output variable discrete_hessians_neigh). - // More precisely, discrete_hessians[i][q_point] stores - // $H_h(\varphi_i)(x_q)$, where $\varphi_i$ is a basis function with support - // on cell, while discrete_hessians_neigh[face_no][i][q_point] stores - // $H_h(\varphi_i)(x_q)$, where $\varphi_i$ is a basis function of the neighboring - // cell adjacent to the face face=cell->face(face_no). - void assemble_local_matrix(const FEValues &fe_values_lift, const unsigned int n_q_points, FullMatrix &local_matrix); - - void compute_discrete_hessians(const typename DoFHandler::active_cell_iterator &cell, - const typename DoFHandler::active_cell_iterator &cell_lift, - std::vector>> &discrete_hessians, - std::vector>>> &discrete_hessians_neigh); - - Triangulation triangulation; + // compute_discrete_hessians computes the required discrete + // Hessians: the discrete Hessians of the basis functions with support on + // the current cell (stored in the output variable + // discrete_hessians) and the basis functions with support on a + // neighbor of the current cell (stored in the output variable + // discrete_hessians_neigh). More precisely, + // discrete_hessians[i][q_point] stores $H_h(\varphi_i)(x_q)$, + // where $\varphi_i$ is a basis function with support on cell, while + // discrete_hessians_neigh[face_no][i][q_point] stores + // $H_h(\varphi_i)(x_q)$, where $\varphi_i$ is a basis function of the + // neighboring cell adjacent to the face + // face=cell->face(face_no). + void assemble_local_matrix(const FEValues &fe_values_lift, + const unsigned int n_q_points, + FullMatrix & local_matrix); + + void compute_discrete_hessians( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator &cell_lift, + std::vector>> & discrete_hessians, + std::vector>>> + &discrete_hessians_neigh); + + Triangulation triangulation; FE_DGQ fe; DoFHandler dof_handler; @@ -131,7 +135,6 @@ namespace Step82 // respectively. double penalty_jump_grad; double penalty_jump_val; - }; @@ -144,33 +147,42 @@ namespace Step82 class RightHandSide : public Function { public: - RightHandSide () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const override; + RightHandSide() + : Function() + {} + virtual double value(const Point & p, + const unsigned int component = 0) const override; }; template - double RightHandSide::value (const Point &p, - const unsigned int /*component*/) const + double RightHandSide::value(const Point &p, + const unsigned int /*component*/) const { double return_value = 0.0; - if (dim==2){ - - return_value = 24.0*std::pow(p(1)*(1.0-p(1)),2)+ - +24.0*std::pow(p(0)*(1.0-p(0)),2) - +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(1)+12.0*p(1)*p(1)); - - } else if (dim==3){ - - return_value = 24.0*std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2) - +24.0*std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2) - +24.0*std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2) - +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(1)+12.0*p(1)*p(1))*std::pow(p(2)*(1.0-p(2)),2) - +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(2)+12.0*p(2)*p(2))*std::pow(p(1)*(1.0-p(1)),2) - +2.0*(2.0-12.0*p(1)+12.0*p(1)*p(1))*(2.0-12.0*p(2)+12.0*p(2)*p(2))*std::pow(p(0)*(1.0-p(0)),2); - - } + if (dim == 2) + { + return_value = 24.0 * std::pow(p(1) * (1.0 - p(1)), 2) + + +24.0 * std::pow(p(0) * (1.0 - p(0)), 2) + + 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * + (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)); + } + else if (dim == 3) + { + return_value = + 24.0 * std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2) + + 24.0 * std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2) + + 24.0 * std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2) + + 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * + (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * + std::pow(p(2) * (1.0 - p(2)), 2) + + 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * + (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * + std::pow(p(1) * (1.0 - p(1)), 2) + + 2.0 * (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * + (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * + std::pow(p(0) * (1.0 - p(0)), 2); + } return return_value; } @@ -178,86 +190,128 @@ namespace Step82 // This class implement the manufactured (exact) solution $u$. To compute the - // errors, we need the value of $u$ as well as its gradient and its Hessian. + // errors, we need the value of $u$ as well as its gradient and its Hessian. template class ExactSolution : public Function { public: - ExactSolution () : Function() {} + ExactSolution() + : Function() + {} - virtual double value (const Point &p, - const unsigned int component = 0) const override; + virtual double value(const Point & p, + const unsigned int component = 0) const override; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const override; - - virtual SymmetricTensor<2,dim> hessian (const Point &p, - const unsigned int component = 0) const override; + virtual Tensor<1, dim> + gradient(const Point & p, + const unsigned int component = 0) const override; + + virtual SymmetricTensor<2, dim> + hessian(const Point & p, + const unsigned int component = 0) const override; }; template - double ExactSolution::value (const Point &p, - const unsigned int /*component*/) const + double ExactSolution::value(const Point &p, + const unsigned int /*component*/) const { double return_value = 0.0; - - if (dim==2){ - return_value = std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2); - } else if (dim==3){ - return_value = std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2); - } - + + if (dim == 2) + { + return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2); + } + else if (dim == 3) + { + return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)) * + p(2) * (1.0 - p(2)), + 2); + } + return return_value; } template - Tensor<1,dim> ExactSolution::gradient (const Point &p, - const unsigned int /*component*/) const + Tensor<1, dim> + ExactSolution::gradient(const Point &p, + const unsigned int /*component*/) const { - Tensor<1,dim> return_gradient; + Tensor<1, dim> return_gradient; return_gradient = 0.0; - if (dim==2){ - return_gradient[0] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * std::pow(p(1)*(1.0-p(1)),2); - return_gradient[1] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * std::pow(p(0)*(1.0-p(0)),2); - } else if (dim==3){ - return_gradient[0] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2); - return_gradient[1] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2); - return_gradient[2] = (2.0*p(2)-6.0*std::pow(p(2),2)+4.0*std::pow(p(2),3)) * std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2); - } - + if (dim == 2) + { + return_gradient[0] = + (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) * + std::pow(p(1) * (1.0 - p(1)), 2); + return_gradient[1] = + (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) * + std::pow(p(0) * (1.0 - p(0)), 2); + } + else if (dim == 3) + { + return_gradient[0] = + (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) * + std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2); + return_gradient[1] = + (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) * + std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2); + return_gradient[2] = + (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) * + std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2); + } + return return_gradient; } template - SymmetricTensor<2,dim> ExactSolution::hessian (const Point &p, - const unsigned int /*component*/) const + SymmetricTensor<2, dim> + ExactSolution::hessian(const Point &p, + const unsigned int /*component*/) const { - SymmetricTensor<2,dim> return_hessian; + SymmetricTensor<2, dim> return_hessian; return_hessian = 0.0; - - if (dim==2){ - return_hessian[0][0] = (2.0-12.0*p(0)+12.0*p(0)*p(0)) * std::pow(p(1)*(1.0-p(1)),2); - return_hessian[0][1] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(1)-6.0*std::pow(p(1),2) - +4.0*std::pow(p(1),3)); - return_hessian[1][1] = (2.0-12.0*p(1)+12.0*p(1)*p(1)) * std::pow(p(0)*(1.0-p(0)),2); - } else if (dim==3){ - return_hessian[0][0] = (2.0-12.0*p(0)+12.0*p(0)*p(0)) * std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2); - return_hessian[0][1] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(1)-6.0*std::pow(p(1),2) - +4.0*std::pow(p(1),3)) * std::pow(p(2)*(1.0-p(2)),2); - return_hessian[0][2] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(2)-6.0*std::pow(p(2),2) - +4.0*std::pow(p(2),3)) * std::pow(p(1)*(1.0-p(1)),2); - return_hessian[1][1] = (2.0-12.0*p(1)+12.0*p(1)*p(1)) * std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2); - return_hessian[1][2] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * (2.0*p(2)-6.0*std::pow(p(2),2) - +4.0*std::pow(p(2),3)) * std::pow(p(0)*(1.0-p(0)),2); - return_hessian[2][2] = (2.0-12.0*p(2)+12.0*p(2)*p(2)) * std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2); - } + + if (dim == 2) + { + return_hessian[0][0] = (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * + std::pow(p(1) * (1.0 - p(1)), 2); + return_hessian[0][1] = + (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) * + (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)); + return_hessian[1][1] = (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * + std::pow(p(0) * (1.0 - p(0)), 2); + } + else if (dim == 3) + { + return_hessian[0][0] = + (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * + std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2); + return_hessian[0][1] = + (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) * + (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) * + std::pow(p(2) * (1.0 - p(2)), 2); + return_hessian[0][2] = + (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) * + (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) * + std::pow(p(1) * (1.0 - p(1)), 2); + return_hessian[1][1] = + (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * + std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2); + return_hessian[1][2] = + (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) * + (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) * + std::pow(p(0) * (1.0 - p(0)), 2); + return_hessian[2][2] = + (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * + std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2); + } return return_hessian; } @@ -272,13 +326,15 @@ namespace Step82 // spaces, we associate the corresponding DoF handlers to the triangulation, // and we set the two penalty coefficients. template - BiLaplacianLDGLift::BiLaplacianLDGLift (const unsigned int fe_degree,double penalty_jump_grad, double penalty_jump_val): - fe(fe_degree), - dof_handler(triangulation), - fe_lift(FE_DGQ(fe_degree),dim*dim), - dof_handler_lift(triangulation), - penalty_jump_grad(penalty_jump_grad), - penalty_jump_val(penalty_jump_val) + BiLaplacianLDGLift::BiLaplacianLDGLift(const unsigned int fe_degree, + double penalty_jump_grad, + double penalty_jump_val) + : fe(fe_degree) + , dof_handler(triangulation) + , fe_lift(FE_DGQ(fe_degree), dim * dim) + , dof_handler_lift(triangulation) + , penalty_jump_grad(penalty_jump_grad) + , penalty_jump_val(penalty_jump_val) {} @@ -286,7 +342,7 @@ namespace Step82 // @sect4{BiLaplacianLDGLift::make_grid} // To build a mesh for $\Omega=(0,1)^d$, we simply call the function - // GridGenerator::hyper_cube and then refine it using + // GridGenerator::hyper_cube and then refine it using // refine_global. The number of refinements is hard-coded // here. template @@ -294,67 +350,76 @@ namespace Step82 { std::cout << "Building the mesh............." << std::endl; - GridGenerator::hyper_cube(triangulation,0.0,1.0); + GridGenerator::hyper_cube(triangulation, 0.0, 1.0); triangulation.refine_global(3); - - std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl; + + std::cout << "Number of active cells: " << triangulation.n_active_cells() + << std::endl; } // @sect4{BiLaplacianLDGLift::setup_system} - // In the following function, we set up the degrees of freedom, the sparsity pattern, - // the size of the matrix $A$, and the size of the solution and right-hand side vectors - // $\boldsymbol{U}$ and $\boldsymbol{F}$. For the sparsity pattern, we cannot directly - // use the function DoFTools::make_flux_sparsity_pattern (as we would do for - // instance for the SIPG method) because we need to take into account the interactions - // of a neighboring cell with another neighboring cell as described in the introduction. - // The extended sparsity pattern is build by iterating over all the active cells. For - // the current cell, we collect all its degrees of freedom as well as the degrees of - // freedom of all its neighboring cells, and then couple everything with everything. + // In the following function, we set up the degrees of freedom, the sparsity + // pattern, the size of the matrix $A$, and the size of the solution and + // right-hand side vectors + // $\boldsymbol{U}$ and $\boldsymbol{F}$. For the sparsity pattern, we cannot + // directly use the function DoFTools::make_flux_sparsity_pattern + // (as we would do for instance for the SIPG method) because we need to take + // into account the interactions of a neighboring cell with another + // neighboring cell as described in the introduction. The extended sparsity + // pattern is build by iterating over all the active cells. For the current + // cell, we collect all its degrees of freedom as well as the degrees of + // freedom of all its neighboring cells, and then couple everything with + // everything. template void BiLaplacianLDGLift::setup_system() { dof_handler.distribute_dofs(fe); dof_handler_lift.distribute_dofs(fe_lift); - - std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; - DynamicSparsityPattern dsp(dof_handler.n_dofs(),dof_handler.n_dofs()); + std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; - const auto dofs_per_cell = fe.dofs_per_cell; + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); - for (const auto cell : dof_handler.active_cell_iterators()){ + const auto dofs_per_cell = fe.dofs_per_cell; - std::vector dofs(dofs_per_cell); - cell->get_dof_indices(dofs); + for (const auto cell : dof_handler.active_cell_iterators()) + { + std::vector dofs(dofs_per_cell); + cell->get_dof_indices(dofs); - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f){ - if (!cell->face(f)->at_boundary()){ - const auto neighbor_cell = cell->neighbor(f); + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (!cell->face(f)->at_boundary()) + { + const auto neighbor_cell = cell->neighbor(f); - std::vector tmp(dofs_per_cell); - neighbor_cell->get_dof_indices(tmp); + std::vector tmp(dofs_per_cell); + neighbor_cell->get_dof_indices(tmp); - dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp)); - } - } + dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp)); + } + } - for (const auto i : dofs){ - for (const auto j : dofs){ - dsp.add(i, j); - dsp.add(j, i); - } + for (const auto i : dofs) + { + for (const auto j : dofs) + { + dsp.add(i, j); + dsp.add(j, i); + } + } } - } sparsity_pattern.copy_from(dsp); std::ofstream out("sparsity_pattern.svg"); sparsity_pattern.print_svg(out); - + matrix.reinit(sparsity_pattern); rhs.reinit(dof_handler.n_dofs()); @@ -374,7 +439,7 @@ namespace Step82 assemble_matrix(); assemble_rhs(); - + std::cout << "Done. " << std::endl; } @@ -390,285 +455,406 @@ namespace Step82 { matrix = 0; - QGauss quad(fe.degree+1); - QGauss quad_face(fe.degree+1); + QGauss quad(fe.degree + 1); + QGauss quad_face(fe.degree + 1); - const unsigned int n_q_points = quad.size(); - const unsigned int n_q_points_face = quad_face.size(); + const unsigned int n_q_points = quad.size(); + const unsigned int n_q_points_face = quad_face.size(); - FEValues fe_values (fe, quad, update_hessians | - update_JxW_values); + FEValues fe_values(fe, quad, update_hessians | update_JxW_values); - FEFaceValues fe_face (fe, quad_face, update_values | - update_gradients | - update_normal_vectors); + FEFaceValues fe_face( + fe, quad_face, update_values | update_gradients | update_normal_vectors); - FEFaceValues fe_face_neighbor (fe, quad_face, update_values | - update_gradients | - update_normal_vectors); + FEFaceValues fe_face_neighbor( + fe, quad_face, update_values | update_gradients | update_normal_vectors); - const unsigned int n_dofs = fe_values.dofs_per_cell; + const unsigned int n_dofs = fe_values.dofs_per_cell; std::vector local_dof_indices(n_dofs), - local_dof_indices_neighbor (n_dofs), - local_dof_indices_neighbor_2 (n_dofs); + local_dof_indices_neighbor(n_dofs), local_dof_indices_neighbor_2(n_dofs); // As indicated in the introduction, the following matrices are used for // the contributions of the products of the discrete Hessians. - FullMatrix stiffness_matrix_cc (n_dofs,n_dofs); // interactions cell / cell - FullMatrix stiffness_matrix_cn (n_dofs,n_dofs); // interactions cell / neighboor - FullMatrix stiffness_matrix_nc (n_dofs,n_dofs); // interactions neighboor / cell - FullMatrix stiffness_matrix_nn (n_dofs,n_dofs); // interactions neighboor / neighboor - FullMatrix stiffness_matrix_n1n2 (n_dofs,n_dofs); // interactions neighboor_1 / neighboor_2 - FullMatrix stiffness_matrix_n2n1 (n_dofs,n_dofs); // interactions neighboor_2 / neighboor_1 + FullMatrix stiffness_matrix_cc(n_dofs, + n_dofs); // interactions cell / cell + FullMatrix stiffness_matrix_cn( + n_dofs, n_dofs); // interactions cell / neighboor + FullMatrix stiffness_matrix_nc( + n_dofs, n_dofs); // interactions neighboor / cell + FullMatrix stiffness_matrix_nn( + n_dofs, n_dofs); // interactions neighboor / neighboor + FullMatrix stiffness_matrix_n1n2( + n_dofs, n_dofs); // interactions neighboor_1 / neighboor_2 + FullMatrix stiffness_matrix_n2n1( + n_dofs, n_dofs); // interactions neighboor_2 / neighboor_1 // The following matrices are used for the contributions of the two // penalty terms. - FullMatrix ip_matrix_cc (n_dofs,n_dofs); // interactions cell / cell - FullMatrix ip_matrix_cn (n_dofs,n_dofs); // interactions cell / neighboor - FullMatrix ip_matrix_nc (n_dofs,n_dofs); // interactions neighboor / cell - FullMatrix ip_matrix_nn (n_dofs,n_dofs); // interactions neighboor / neighboor - - std::vector>> discrete_hessians (n_dofs, std::vector>(n_q_points) ); - std::vector>>> discrete_hessians_neigh (GeometryInfo::faces_per_cell, discrete_hessians); - - Tensor<2,dim> H_i,H_j; - Tensor<2,dim> H_i_neigh,H_j_neigh; - Tensor<2,dim> H_i_neigh2,H_j_neigh2; - - double mesh_inv,mesh3_inv; - bool at_boundary,at_boundary_2; + FullMatrix ip_matrix_cc(n_dofs, n_dofs); // interactions cell / cell + FullMatrix ip_matrix_cn(n_dofs, + n_dofs); // interactions cell / neighboor + FullMatrix ip_matrix_nc(n_dofs, + n_dofs); // interactions neighboor / cell + FullMatrix ip_matrix_nn( + n_dofs, n_dofs); // interactions neighboor / neighboor + + std::vector>> discrete_hessians( + n_dofs, std::vector>(n_q_points)); + std::vector>>> + discrete_hessians_neigh(GeometryInfo::faces_per_cell, + discrete_hessians); + + Tensor<2, dim> H_i, H_j; + Tensor<2, dim> H_i_neigh, H_j_neigh; + Tensor<2, dim> H_i_neigh2, H_j_neigh2; + + double mesh_inv, mesh3_inv; + bool at_boundary, at_boundary_2; unsigned int face_no_neighbor = 0; - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - typename DoFHandler::active_cell_iterator neighbor_cell,neighbor_cell_2; - - typename DoFHandler::active_cell_iterator cell_lift = dof_handler_lift.begin_active(); - - for (; cell != endc; ++cell, ++cell_lift){ - - fe_values.reinit(cell); - cell->get_dof_indices (local_dof_indices); - - // We now compute all the discrete Hessians that are not vanishing - // on the current cell, i.e., the discrete Hessian of all the basis - // functions with support on the current cell or on one of its neighbors. - compute_discrete_hessians(cell,cell_lift, - discrete_hessians,discrete_hessians_neigh); - - // First, we compute and add the interactions of the degrees of freedom - // of the current cell. - stiffness_matrix_cc = 0; - for (unsigned int q=0; q::faces_per_cell; ++face_no){ - const typename DoFHandler::face_iterator face=cell->face(face_no); - - at_boundary = face->at_boundary(); - - if (!at_boundary){ // nothing to be done if boundary face (the liftings of the Dirichlet BCs are accounted for in the assembly of the RHS; in fact, nothing to be done in this program since we prescribe homogeneous BCs) - - neighbor_cell =cell->neighbor(face_no); - neighbor_cell->get_dof_indices (local_dof_indices_neighbor); - - stiffness_matrix_cn=0; - stiffness_matrix_nc=0; - stiffness_matrix_nn=0; - for (unsigned int q=0; q::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + + typename DoFHandler::active_cell_iterator neighbor_cell, + neighbor_cell_2; + + typename DoFHandler::active_cell_iterator cell_lift = + dof_handler_lift.begin_active(); + + for (; cell != endc; ++cell, ++cell_lift) + { + fe_values.reinit(cell); + cell->get_dof_indices(local_dof_indices); + + // We now compute all the discrete Hessians that are not vanishing + // on the current cell, i.e., the discrete Hessian of all the basis + // functions with support on the current cell or on one of its + // neighbors. + compute_discrete_hessians(cell, + cell_lift, + discrete_hessians, + discrete_hessians_neigh); + + // First, we compute and add the interactions of the degrees of freedom + // of the current cell. + stiffness_matrix_cc = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { const double dx = fe_values.JxW(q); - - for (unsigned int i=0; i::faces_per_cell-1; ++face_no){ - const typename DoFHandler::face_iterator face=cell->face(face_no); - - at_boundary = face->at_boundary(); - - if (!at_boundary){ // nothing to be done if boundary face (the liftings of the Dirichlet BCs are accounted for in the assembly of the RHS; in fact, nothing to be done in this program since we prescribe homogeneous BCs) - - for (unsigned int face_no_2=face_no+1; face_no_2 < GeometryInfo::faces_per_cell; ++face_no_2){ - const typename DoFHandler::face_iterator face_2=cell->face(face_no_2); - at_boundary_2 = face_2->at_boundary(); + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + H_i = discrete_hessians[i][q]; + H_j = discrete_hessians[j][q]; - if (!at_boundary_2){ - - neighbor_cell = cell->neighbor(face_no); - neighbor_cell->get_dof_indices (local_dof_indices_neighbor); - neighbor_cell_2 = cell->neighbor(face_no_2); - neighbor_cell_2->get_dof_indices (local_dof_indices_neighbor_2); - - stiffness_matrix_n1n2=0; - stiffness_matrix_n2n1=0; - - for (unsigned int q=0; q::faces_per_cell; ++face_no){ - const typename DoFHandler::face_iterator face=cell->face(face_no); - - mesh_inv = 1.0/face->diameter(); // h_e^{-1} - mesh3_inv = 1.0/std::pow(face->diameter(),3); // ĥ_e^{-3} - - fe_face.reinit(cell,face_no); - - ip_matrix_cc = 0; // filled in any case (boundary or interior face) - - at_boundary = face->at_boundary(); - if (at_boundary){ - - for (unsigned int q=0; qneighbor(face_no); - face_no_neighbor = cell->neighbor_of_neighbor (face_no); - - if(neighbor_cell->id().operator<(cell->id())){ //we need to have a global way to compare the cells in order to not calculate the same jump term twice - continue; // skip this face (already considered) - } else{ - - fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor); - neighbor_cell->get_dof_indices (local_dof_indices_neighbor); - - ip_matrix_cn = 0; - ip_matrix_nc = 0; - ip_matrix_nn = 0; - - for (unsigned int q=0; q::faces_per_cell; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + + at_boundary = face->at_boundary(); + + if (!at_boundary) + { // nothing to be done if boundary face (the liftings of the + // Dirichlet BCs are accounted for in the assembly of the RHS; + // in fact, nothing to be done in this program since we + // prescribe homogeneous BCs) + + neighbor_cell = cell->neighbor(face_no); + neighbor_cell->get_dof_indices(local_dof_indices_neighbor); + + stiffness_matrix_cn = 0; + stiffness_matrix_nc = 0; + stiffness_matrix_nn = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double dx = fe_values.JxW(q); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + H_i = discrete_hessians[i][q]; + H_j = discrete_hessians[j][q]; + + H_i_neigh = discrete_hessians_neigh[face_no][i][q]; + H_j_neigh = discrete_hessians_neigh[face_no][j][q]; + + stiffness_matrix_cn(i, j) += + dx * scalar_product(H_j_neigh, H_i); + stiffness_matrix_nc(i, j) += + dx * scalar_product(H_j, H_i_neigh); + stiffness_matrix_nn(i, j) += + dx * scalar_product(H_j_neigh, H_i_neigh); + } + } + } + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + matrix(local_dof_indices[i], + local_dof_indices_neighbor[j]) += + stiffness_matrix_cn(i, j); + matrix(local_dof_indices_neighbor[i], + local_dof_indices[j]) += + stiffness_matrix_nc(i, j); + matrix(local_dof_indices_neighbor[i], + local_dof_indices_neighbor[j]) += + stiffness_matrix_nn(i, j); + } + } + + } // boundary check + } // for face + + // We now compute and add the interactions of the degrees of freedom of + // a neighboring cells with those of another neighboring cell (this is + // where we need the extended sparsity pattern). + for (unsigned int face_no = 0; + face_no < GeometryInfo::faces_per_cell - 1; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + + at_boundary = face->at_boundary(); + + if (!at_boundary) + { // nothing to be done if boundary face (the liftings of the + // Dirichlet BCs are accounted for in the assembly of the RHS; + // in fact, nothing to be done in this program since we + // prescribe homogeneous BCs) + + + for (unsigned int face_no_2 = face_no + 1; + face_no_2 < GeometryInfo::faces_per_cell; + ++face_no_2) + { + const typename DoFHandler::face_iterator face_2 = + cell->face(face_no_2); + at_boundary_2 = face_2->at_boundary(); + + if (!at_boundary_2) + { + neighbor_cell = cell->neighbor(face_no); + neighbor_cell->get_dof_indices( + local_dof_indices_neighbor); + neighbor_cell_2 = cell->neighbor(face_no_2); + neighbor_cell_2->get_dof_indices( + local_dof_indices_neighbor_2); + + stiffness_matrix_n1n2 = 0; + stiffness_matrix_n2n1 = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double dx = fe_values.JxW(q); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + H_i_neigh = + discrete_hessians_neigh[face_no][i][q]; + H_j_neigh = + discrete_hessians_neigh[face_no][j][q]; + + H_i_neigh2 = + discrete_hessians_neigh[face_no_2][i][q]; + H_j_neigh2 = + discrete_hessians_neigh[face_no_2][j][q]; + + stiffness_matrix_n1n2(i, j) += + dx * + scalar_product(H_j_neigh2, H_i_neigh); + stiffness_matrix_n2n1(i, j) += + dx * + scalar_product(H_j_neigh, H_i_neigh2); + } + } + } + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + matrix(local_dof_indices_neighbor[i], + local_dof_indices_neighbor_2[j]) += + stiffness_matrix_n1n2(i, j); + matrix(local_dof_indices_neighbor_2[i], + local_dof_indices_neighbor[j]) += + stiffness_matrix_n2n1(i, j); + } + } + } // boundary check face_2 + } // for face_2 + } // boundary check face_1 + } // for face_1 + + + // Finally, we compute and add the two penalty terms. + for (unsigned int face_no = 0; + face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + + mesh_inv = 1.0 / face->diameter(); // h_e^{-1} + mesh3_inv = 1.0 / std::pow(face->diameter(), 3); // ĥ_e^{-3} + + fe_face.reinit(cell, face_no); + + ip_matrix_cc = 0; // filled in any case (boundary or interior face) + + at_boundary = face->at_boundary(); + if (at_boundary) + { + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face.JxW(q); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + ip_matrix_cc(i, j) += penalty_jump_grad * mesh_inv * + dx * + fe_face.shape_grad(j, q) * + fe_face.shape_grad(i, q); + ip_matrix_cc(i, j) += penalty_jump_val * mesh3_inv * + dx * + fe_face.shape_value(j, q) * + fe_face.shape_value(i, q); + } + } + } + } + else + { // interior face + + neighbor_cell = cell->neighbor(face_no); + face_no_neighbor = cell->neighbor_of_neighbor(face_no); + + if (neighbor_cell->id().operator<(cell->id())) + { // we need to have a global way to compare the cells in + // order to not calculate the same jump term twice + continue; // skip this face (already considered) + } + else + { + fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor); + neighbor_cell->get_dof_indices(local_dof_indices_neighbor); + + ip_matrix_cn = 0; + ip_matrix_nc = 0; + ip_matrix_nn = 0; + + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face.JxW(q); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + ip_matrix_cc(i, j) += penalty_jump_grad * + mesh_inv * dx * + fe_face.shape_grad(j, q) * + fe_face.shape_grad(i, q); + ip_matrix_cc(i, j) += + penalty_jump_val * mesh3_inv * dx * + fe_face.shape_value(j, q) * + fe_face.shape_value(i, q); + + ip_matrix_cn(i, j) -= + penalty_jump_grad * mesh_inv * dx * + fe_face_neighbor.shape_grad(j, q) * + fe_face.shape_grad(i, q); + ip_matrix_cn(i, j) -= + penalty_jump_val * mesh3_inv * dx * + fe_face_neighbor.shape_value(j, q) * + fe_face.shape_value(i, q); + + ip_matrix_nc(i, j) -= + penalty_jump_grad * mesh_inv * dx * + fe_face.shape_grad(j, q) * + fe_face_neighbor.shape_grad(i, q); + ip_matrix_nc(i, j) -= + penalty_jump_val * mesh3_inv * dx * + fe_face.shape_value(j, q) * + fe_face_neighbor.shape_value(i, q); + + ip_matrix_nn(i, j) += + penalty_jump_grad * mesh_inv * dx * + fe_face_neighbor.shape_grad(j, q) * + fe_face_neighbor.shape_grad(i, q); + ip_matrix_nn(i, j) += + penalty_jump_val * mesh3_inv * dx * + fe_face_neighbor.shape_value(j, q) * + fe_face_neighbor.shape_value(i, q); + } + } + } + } // face not visited yet + + } // boundary check + + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + matrix(local_dof_indices[i], local_dof_indices[j]) += + ip_matrix_cc(i, j); + } + } - for (unsigned int i=0; i::assemble_rhs() { rhs = 0; - - QGauss quad(fe.degree+1); - FEValues fe_values(fe, quad, update_values | - update_quadrature_points | - update_JxW_values); - const unsigned int n_dofs = fe_values.dofs_per_cell; + QGauss quad(fe.degree + 1); + FEValues fe_values( + fe, quad, update_values | update_quadrature_points | update_JxW_values); + + const unsigned int n_dofs = fe_values.dofs_per_cell; const unsigned int n_quad_pts = quad.size(); const RightHandSide right_hand_side; - Vector local_rhs(n_dofs); + Vector local_rhs(n_dofs); std::vector local_dof_indices(n_dofs); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell!=endc; ++cell){ + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); - fe_values.reinit(cell); - cell->get_dof_indices(local_dof_indices); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + cell->get_dof_indices(local_dof_indices); - local_rhs = 0; - for (unsigned int q=0; q void BiLaplacianLDGLift::compute_errors() - { - - double error_H2 = 0; // sqrt( ||D_h^2(u-u_h)||_{L^2(Omega)}^2 + ||h^{-1/2}[grad_h(u-u_h)]||_{L^2(Sigma)}^2 + ||h^{-3/2}[u-u_h]||_{L^2(Sigma)}^2 ) - double error_H1 = 0; // sqrt( ||grad_h(u-u_h)||_{L^2(Omega)}^2 + ||h^{-1/2}[u-u_h]||_{L^2(Sigma)}^2 ) + { + double error_H2 = 0; // sqrt( ||D_h^2(u-u_h)||_{L^2(Omega)}^2 + + // ||h^{-1/2}[grad_h(u-u_h)]||_{L^2(Sigma)}^2 + + // ||h^{-3/2}[u-u_h]||_{L^2(Sigma)}^2 ) + double error_H1 = 0; // sqrt( ||grad_h(u-u_h)||_{L^2(Omega)}^2 + + // ||h^{-1/2}[u-u_h]||_{L^2(Sigma)}^2 ) double error_L2 = 0; // ||u-u_h||_{L^2(Omega)} - - QGauss quad(fe.degree+1); - QGauss quad_face(fe.degree+1); - - FEValues fe_values (fe, quad, update_values | - update_gradients | - update_hessians | - update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face (fe, quad_face, update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face_neighbor (fe, quad_face, update_values | - update_gradients); - - const unsigned int n_q_points = quad.size(); - const unsigned int n_q_points_face = quad_face.size(); - + + QGauss quad(fe.degree + 1); + QGauss quad_face(fe.degree + 1); + + FEValues fe_values(fe, + quad, + update_values | update_gradients | update_hessians | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face(fe, + quad_face, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_neighbor(fe, + quad_face, + update_values | update_gradients); + + const unsigned int n_q_points = quad.size(); + const unsigned int n_q_points_face = quad_face.size(); + // We introduce some variables for the exact solution const ExactSolution u_exact; - double u_exact_q; - Tensor<1,dim> u_exact_grad_q; - + double u_exact_q; + Tensor<1, dim> u_exact_grad_q; + // and for the approximate solution - std::vector solution_values_cell(n_q_points); - std::vector> solution_gradients_cell(n_q_points); - std::vector> solution_hessians_cell(n_q_points); - - std::vector solution_values(n_q_points_face); - std::vector solution_values_neigh(n_q_points_face); - std::vector> solution_gradients(n_q_points_face); - std::vector> solution_gradients_neigh(n_q_points_face); - + std::vector solution_values_cell(n_q_points); + std::vector> solution_gradients_cell(n_q_points); + std::vector> solution_hessians_cell(n_q_points); + + std::vector solution_values(n_q_points_face); + std::vector solution_values_neigh(n_q_points_face); + std::vector> solution_gradients(n_q_points_face); + std::vector> solution_gradients_neigh(n_q_points_face); + double mesh_inv; double mesh3_inv; - bool at_boundary; - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - + bool at_boundary; + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator neighbor_cell; - unsigned int face_no_neighbor = 0; - - for (; cell!=endc; ++cell){ - - fe_values.reinit (cell); - - fe_values.get_function_values(solution,solution_values_cell); - fe_values.get_function_gradients(solution,solution_gradients_cell); - fe_values.get_function_hessians(solution,solution_hessians_cell); - - // We first add the bulk terms. - for (unsigned int q=0; q::faces_per_cell; ++face_no){ - const typename DoFHandler::face_iterator face=cell->face(face_no); - - mesh_inv = 1.0/face->diameter(); // h^{-1} - mesh3_inv = 1.0/std::pow(face->diameter(),3); // h^{-3} - - fe_face.reinit(cell,face_no); - - fe_face.get_function_values(solution,solution_values); - fe_face.get_function_gradients(solution,solution_gradients); - - at_boundary = face->at_boundary(); - if (at_boundary){ - - for (unsigned int q=0; qneighbor(face_no); - face_no_neighbor = cell->neighbor_of_neighbor (face_no); - - if(neighbor_cell->id().operator<(cell->id())){ // we need to have a global way to compare the cells in order to not calculate the same jump term twice - continue; // skip this face (already considered) - } else{ - - fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor); - - fe_face.get_function_values(solution,solution_values); - fe_face_neighbor.get_function_values(solution,solution_values_neigh); - fe_face.get_function_gradients(solution,solution_gradients); - fe_face_neighbor.get_function_gradients(solution,solution_gradients_neigh); - - for (unsigned int q=0; qbulk terms. + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double dx = fe_values.JxW(q); + + error_H2 += dx * (u_exact.hessian(fe_values.quadrature_point(q)) - + solution_hessians_cell[q]) + .norm_square(); + error_H1 += dx * (u_exact.gradient(fe_values.quadrature_point(q)) - + solution_gradients_cell[q]) + .norm_square(); + error_L2 += + dx * std::pow(u_exact.value(fe_values.quadrature_point(q)) - + solution_values_cell[q], + 2); + } // for quad + + // We then add the face contributions. + for (unsigned int face_no = 0; + face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + + mesh_inv = 1.0 / face->diameter(); // h^{-1} + mesh3_inv = 1.0 / std::pow(face->diameter(), 3); // h^{-3} + + fe_face.reinit(cell, face_no); + + fe_face.get_function_values(solution, solution_values); + fe_face.get_function_gradients(solution, solution_gradients); + + at_boundary = face->at_boundary(); + if (at_boundary) + { + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face.JxW(q); + u_exact_q = u_exact.value(fe_face.quadrature_point(q)); + u_exact_grad_q = + u_exact.gradient(fe_face.quadrature_point(q)); + + error_H2 += + dx * mesh_inv * + (u_exact_grad_q - solution_gradients[q]).norm_square(); + error_H2 += dx * mesh3_inv * + std::pow(u_exact_q - solution_values[q], 2); + error_H1 += dx * mesh_inv * + std::pow(u_exact_q - solution_values[q], 2); + } + } + else + { // interior face + + neighbor_cell = cell->neighbor(face_no); + face_no_neighbor = cell->neighbor_of_neighbor(face_no); + + if (neighbor_cell->id().operator<(cell->id())) + { // we need to have a global way to compare the cells in + // order to not calculate the same jump term twice + continue; // skip this face (already considered) + } + else + { + fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor); + + fe_face.get_function_values(solution, solution_values); + fe_face_neighbor.get_function_values(solution, + solution_values_neigh); + fe_face.get_function_gradients(solution, + solution_gradients); + fe_face_neighbor.get_function_gradients( + solution, solution_gradients_neigh); + + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face.JxW(q); + + // To compute the jump term, we use the fact that + // $\jump{u}=0$ and + // $\jump{\nabla u}=\mathbf{0}$ since $u\in + // H^2(\Omega)$. + error_H2 += + dx * mesh_inv * + (solution_gradients_neigh[q] - solution_gradients[q]) + .norm_square(); + error_H2 += dx * mesh3_inv * + std::pow(solution_values_neigh[q] - + solution_values[q], + 2); + error_H1 += dx * mesh_inv * + std::pow(solution_values_neigh[q] - + solution_values[q], + 2); + } + } // face not visited yet + + } // boundary check + + } // for face + + } // for cell + error_H2 = std::sqrt(error_H2); error_H1 = std::sqrt(error_H1); error_L2 = std::sqrt(error_L2); std::cout << "DG H2 norm of the error: " << error_H2 << std::endl; std::cout << "DG H1 norm of the error: " << error_H1 << std::endl; - std::cout << " L2 norm of the error: " << error_L2 << std::endl; + std::cout << " L2 norm of the error: " << error_L2 << std::endl; } @@ -894,11 +1123,11 @@ namespace Step82 { DataOut data_out; data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "solution"); + data_out.add_data_vector(solution, "solution"); data_out.build_patches(); - std::ofstream output ("solution.vtk"); - data_out.write_vtk (output); + std::ofstream output("solution.vtk"); + data_out.write_vtk(output); } @@ -910,23 +1139,30 @@ namespace Step82 // lifting terms. We reiterate that only the basis functions with // support on the current cell are accounting for. template - void BiLaplacianLDGLift::assemble_local_matrix(const FEValues &fe_values_lift, const unsigned int n_q_points, FullMatrix &local_matrix) + void BiLaplacianLDGLift::assemble_local_matrix( + const FEValues &fe_values_lift, + const unsigned int n_q_points, + FullMatrix & local_matrix) { const FEValuesExtractors::Tensor<2> tau_ext(0); const unsigned int n_dofs = fe_values_lift.dofs_per_cell; local_matrix = 0; - for (unsigned int q=0; q - void BiLaplacianLDGLift::compute_discrete_hessians(const typename DoFHandler::active_cell_iterator &cell, - const typename DoFHandler::active_cell_iterator &cell_lift, - std::vector>> &discrete_hessians, - std::vector>>> &discrete_hessians_neigh) + void BiLaplacianLDGLift::compute_discrete_hessians( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator &cell_lift, + std::vector>> & discrete_hessians, + std::vector>>> + &discrete_hessians_neigh) { - QGauss quad(fe.degree+1); - QGauss quad_face(fe.degree+1); + QGauss quad(fe.degree + 1); + QGauss quad_face(fe.degree + 1); - const unsigned int n_q_points = quad.size(); - const unsigned int n_q_points_face = quad_face.size(); + const unsigned int n_q_points = quad.size(); + const unsigned int n_q_points_face = quad_face.size(); // The information we need from the basis functions of // $\mathbb{V}_h$: fe_values is needed to add @@ -957,21 +1196,18 @@ namespace Step82 // fe_face and fe_face_neighbor // are used to compute the right-hand sides for the local // problems. - FEValues fe_values (fe, quad, update_hessians | - update_JxW_values); + FEValues fe_values(fe, quad, update_hessians | update_JxW_values); - FEFaceValues fe_face (fe, quad_face, update_values | - update_gradients | - update_normal_vectors); + FEFaceValues fe_face( + fe, quad_face, update_values | update_gradients | update_normal_vectors); - FEFaceValues fe_face_neighbor (fe, quad_face, update_values | - update_gradients | - update_normal_vectors); + FEFaceValues fe_face_neighbor( + fe, quad_face, update_values | update_gradients | update_normal_vectors); const unsigned int n_dofs = fe_values.dofs_per_cell; - typename DoFHandler<2,dim>::active_cell_iterator neighbor_cell; - unsigned int face_no_neighbor = 0; + typename DoFHandler<2, dim>::active_cell_iterator neighbor_cell; + unsigned int face_no_neighbor = 0; // The information needed from the basis functions // of the finite element space for the lifting terms: @@ -979,186 +1215,246 @@ namespace Step82 // mass matrix (see $\boldsymbol{M}_c$ in the introduction), // while fe_face_lift is used to compute the // right-hand sides (see $\boldsymbol{G}_c$ for $b_e$). - FEValues fe_values_lift (fe_lift, quad, update_values | - update_JxW_values); + FEValues fe_values_lift(fe_lift, + quad, + update_values | update_JxW_values); - FEFaceValues fe_face_lift (fe_lift, quad_face, update_values | - update_gradients | - update_JxW_values); + FEFaceValues fe_face_lift( + fe_lift, quad_face, update_values | update_gradients | update_JxW_values); const FEValuesExtractors::Tensor<2> tau_ext(0); - + const unsigned int n_dofs_lift = fe_values_lift.dofs_per_cell; - FullMatrix local_matrix_lift (n_dofs_lift,n_dofs_lift); + FullMatrix local_matrix_lift(n_dofs_lift, n_dofs_lift); Vector local_rhs_re(n_dofs_lift), local_rhs_be(n_dofs_lift), - coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift), - coeffs_tmp(n_dofs_lift); - + coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift), coeffs_tmp(n_dofs_lift); + SolverControl solver_control(1000, 1e-12); - SolverCG<> solver(solver_control); - - bool at_boundary; + SolverCG<> solver(solver_control); + + bool at_boundary; double factor_avg; // 0.5 for interior faces, 1.0 for boundary faces - + fe_values.reinit(cell); - fe_values_lift.reinit(cell_lift); - + fe_values_lift.reinit(cell_lift); + // We start by assembling the (local) mass matrix used for the computation // of the lifting terms $r_e$ and $b_e$. - assemble_local_matrix(fe_values_lift,n_q_points,local_matrix_lift); - - for (unsigned int i=0; i::faces_per_cell; + ++face_no) + { + discrete_hessians_neigh[face_no][i][q] = 0; + } + } + } - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no){ - discrete_hessians_neigh[face_no][i][q]=0; - } + // In this loop, we compute the discrete Hessian at each quadrature point + // $x_q$ of cell for each basis function supported on + // cell, namely we fill-in the variable + // discrete_hessians[i][q]. For the lifting terms, we need to + // add the contribution of all the faces of cell. + for (unsigned int i = 0; i < n_dofs; ++i) + { + coeffs_re = 0; + coeffs_be = 0; + + for (unsigned int face_no = 0; + face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + + at_boundary = face->at_boundary(); + + // Recall that by convention, the average of a function accross a + // boundary face $e$ reduces to the trace of the function on the + // only element adjacent to $e$, namely there is no factor + // $\frac{1}{2}$. We distinguish between the two cases (the current + // face lies in the interior or on the boundary of the domain) using + // the variable factor_avg. + factor_avg = 0.5; + if (at_boundary) + { + factor_avg = 1.0; + } - } - } + fe_face.reinit(cell, face_no); + fe_face_lift.reinit(cell_lift, face_no); + + local_rhs_re = 0; + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face_lift.JxW(q); + const Tensor<1, dim> normal = fe_face.normal_vector( + q); // same as fe_face_lift.normal_vector(q) + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + local_rhs_re(m) += + factor_avg * dx * + (fe_face_lift[tau_ext].value(m, q) * normal) * + fe_face.shape_grad(i, q); + } + } - // In this loop, we compute the discrete Hessian at each quadrature point $x_q$ - // of cell for each basis function supported on cell, - // namely we fill-in the variable discrete_hessians[i][q]. - // For the lifting terms, we need to add the contribution of all the faces of - // cell. - for (unsigned int i=0; ilocal_rhs_be(m) corresponds to $G_m$ + // introduced in the comments about the implementation of the + // lifting $b_e$ in the case + // $\varphi=\varphi^c$. + local_rhs_be = 0; + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face_lift.JxW(q); + const Tensor<1, dim> normal = fe_face.normal_vector( + q); // same as fe_face_lift.normal_vector(q) + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + local_rhs_be(m) += factor_avg * dx * + fe_face_lift[tau_ext].divergence(m, q) * + normal * fe_face.shape_value(i, q); + } + } - coeffs_re=0; coeffs_be=0; + coeffs_tmp = 0; + solver.solve(local_matrix_lift, + coeffs_tmp, + local_rhs_re, + PreconditionIdentity()); + coeffs_re += coeffs_tmp; + + coeffs_tmp = 0; + solver.solve(local_matrix_lift, + coeffs_tmp, + local_rhs_be, + PreconditionIdentity()); + coeffs_be += coeffs_tmp; + + } // for face + + for (unsigned int q = 0; q < n_q_points; ++q) + { + discrete_hessians[i][q] += fe_values.shape_hessian(i, q); + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + discrete_hessians[i][q] -= + coeffs_re[m] * fe_values_lift[tau_ext].value(m, q); + } - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no){ - const typename DoFHandler::face_iterator face=cell->face(face_no); - - at_boundary = face->at_boundary(); - - // Recall that by convention, the average of a function accross a boundary - // face $e$ reduces to the trace of the function on the only element - // adjacent to $e$, namely there is no factor $\frac{1}{2}$. We distinguish - // between the two cases (the current face lies in the interior or on - // the boundary of the domain) using the variable factor_avg. - factor_avg = 0.5; - if (at_boundary){ - factor_avg = 1.0; - } - - fe_face.reinit(cell,face_no); - fe_face_lift.reinit(cell_lift,face_no); - - local_rhs_re=0; - for (unsigned int q=0; q normal = fe_face.normal_vector(q); // same as fe_face_lift.normal_vector(q) - - for (unsigned int m=0; mlocal_rhs_be(m) corresponds to $G_m$ introduced in the - // comments about the implementation of the lifting $b_e$ in the case - // $\varphi=\varphi^c$. - local_rhs_be=0; - for (unsigned int q=0; q normal = fe_face.normal_vector(q); // same as fe_face_lift.normal_vector(q) - - for (unsigned int m=0; mcell for each basis function supported on a neighboring - // neighbor_cell of cell, namely we fill-in the - // variable discrete_hessians_neigh[face_no][i][q]. + // $x_q$ of cell for each basis function supported on a + // neighboring neighbor_cell of cell, namely we + // fill-in the variable discrete_hessians_neigh[face_no][i][q]. // For the lifting terms, we only need to add the contribution of the - // face adjecent to cell and neighbor_cell. - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no){ - - const typename DoFHandler::face_iterator face=cell->face(face_no); - - at_boundary = face->at_boundary(); - - if (!at_boundary){ // for non-homogeneous Dirichlet BCs, we would need to compute the lifting of the prescribed BC (see Section Possible Extensions for more details) - - neighbor_cell =cell->neighbor(face_no); - face_no_neighbor = cell->neighbor_of_neighbor (face_no); - fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor); - - for (unsigned int i=0; i normal = fe_face_neighbor.normal_vector(q); - - for (unsigned int m=0; mlocal_rhs_be(m) corresponds to $G_m$ introduced in - // the comments about the implementation of the lifting $b_e$ in the case - // $\varphi=\varphi^n$. - local_rhs_be=0; - for (unsigned int q=0; q normal = fe_face_neighbor.normal_vector(q); - - for (unsigned int m=0; mcell and neighbor_cell. + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); - for (unsigned int m=0; mat_boundary(); - } // for dof i - } // boundary check - } // for face + if (!at_boundary) + { // for non-homogeneous Dirichlet BCs, we would need to compute the + // lifting of the prescribed BC (see Section Possible Extensions for + // more details) + + neighbor_cell = cell->neighbor(face_no); + face_no_neighbor = cell->neighbor_of_neighbor(face_no); + fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + coeffs_re = 0; + coeffs_be = 0; + + fe_face_lift.reinit(cell_lift, face_no); + + local_rhs_re = 0; + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face_lift.JxW(q); + const Tensor<1, dim> normal = + fe_face_neighbor.normal_vector(q); + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + local_rhs_re(m) += + 0.5 * dx * + (fe_face_lift[tau_ext].value(m, q) * normal) * + fe_face_neighbor.shape_grad(i, q); + } + } + + // Here, local_rhs_be(m) corresponds to $G_m$ + // introduced in the comments about the implementation of the + // lifting $b_e$ in the case + // $\varphi=\varphi^n$. + local_rhs_be = 0; + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const double dx = fe_face_lift.JxW(q); + const Tensor<1, dim> normal = + fe_face_neighbor.normal_vector(q); + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + local_rhs_be(m) += + 0.5 * dx * fe_face_lift[tau_ext].divergence(m, q) * + normal * fe_face_neighbor.shape_value(i, q); + } + } + + solver.solve(local_matrix_lift, + coeffs_re, + local_rhs_re, + PreconditionIdentity()); + solver.solve(local_matrix_lift, + coeffs_be, + local_rhs_be, + PreconditionIdentity()); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + discrete_hessians_neigh[face_no][i][q] -= + coeffs_re[m] * fe_values_lift[tau_ext].value(m, q); + } + + for (unsigned int m = 0; m < n_dofs_lift; ++m) + { + discrete_hessians_neigh[face_no][i][q] += + coeffs_be[m] * fe_values_lift[tau_ext].value(m, q); + } + } + + } // for dof i + } // boundary check + } // for face } @@ -1168,7 +1464,7 @@ namespace Step82 void BiLaplacianLDGLift::run() { make_grid(); - + setup_system(); assemble_system(); @@ -1178,7 +1474,7 @@ namespace Step82 output_results(); } -} // namespace Step82 +} // namespace Step82 @@ -1187,17 +1483,48 @@ namespace Step82 // The is the main function. We define here the polynomial degree // for the two finite element spaces (for the solution and the two liftings) and // the two penalty coefficients. We can also change the dimension to run the -// code in 3D. +// code in 3D. int main() { - int degree=2; // FE degree for u_h and the two lifting terms + try + { + const unsigned int int degree = + 2; // FE degree for u_h and the two lifting terms + + const double penalty_grad = + 1.0; // penalty coefficient for the jump of the gradients + const double penalty_val = + 1.0; // penalty coefficient for the jump of the values - double penalty_grad=1.0; // penalty coefficient for the jump of the gradients - double penalty_val=1.0; // penalty coefficient for the jump of the values + Step82::BiLaplacianLDGLift<2> problem(degree, penalty_grad, penalty_val); + + problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } - Step82::BiLaplacianLDGLift<2> problem(degree,penalty_grad,penalty_val); - - problem.run(); - return 0; } -- 2.39.5