From e3109ce2d1e75db170b255b736b044e712830514 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 14 Apr 2019 18:08:40 -0600 Subject: [PATCH] Enclose everything in step-61 into its own namespace. --- examples/step-61/step-61.cc | 1600 ++++++++++++++++++----------------- 1 file changed, 809 insertions(+), 791 deletions(-) diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc index d62b125a7b..331876650f 100644 --- a/examples/step-61/step-61.cc +++ b/examples/step-61/step-61.cc @@ -13,12 +13,13 @@ * * --------------------------------------------------------------------- - * Author: Zhuoran Wang + * Author: Zhuoran Wang, Colorado State University, 2018 */ // @sect3{Include files} // This program is based on step-7, step-20 and step-51, -// we add these include files. +// so most of the following header files are familiar. We +// need the following: #include #include #include @@ -58,854 +59,870 @@ #include #include -using namespace dealii; -// @sect3{The WGDarcyEquation class template} - -// We will solve for the numerical pressure in the interior and on faces and -// calculate its $L_2$ error of pressure. In the post-processing step, we will -// calculate $L_2$-errors of velocity and flux. -template -class WGDarcyEquation -{ -public: - WGDarcyEquation(); - void run(); - -private: - void make_grid(); - void setup_system(); - void assemble_system(); - void solve(); - void postprocess(); - void process_solution(); - void output_results() const; - - Triangulation triangulation; - - AffineConstraints constraints; - - FE_RaviartThomas fe_rt; - DoFHandler dof_handler_rt; - - // The finite element system is used for interior and face solutions. - FESystem fe; - DoFHandler dof_handler; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; -}; - -// @sect3{Right hand side, boundary values, and exact solution} - -// Next, we define the coefficient matrix $\mathbf{K}$, -// Dirichlet boundary conditions, the right-hand side $f = 2\pi^2 \sin(\pi x) -// \sin(\pi y)$, and the reference solutions $p = \sin(\pi x) \sin(\pi y) $. -// -// The coefficient matrix $\mathbf{K}$ is the identity matrix as a test example. -template -class Coefficient : public TensorFunction<2, dim> -{ -public: - Coefficient() - : TensorFunction<2, dim>() - {} - - virtual void value_list(const std::vector> &points, - std::vector> & values) const override; -}; - -template -void Coefficient::value_list(const std::vector> &points, - std::vector> & values) const -{ - Assert(points.size() == values.size(), - ExcDimensionMismatch(points.size(), values.size())); - for (unsigned int p = 0; p < points.size(); ++p) - { - values[p].clear(); - for (unsigned int d = 0; d < dim; ++d) - values[p][d][d] = 1; - } -} - -template -class BoundaryValues : public Function -{ -public: - BoundaryValues() - : Function(2) - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; -}; - -template -double BoundaryValues::value(const Point & /*p*/, - const unsigned int /*component*/) const +// Our first step, as always, is to put everything related to this tutorial +// program into its own namespace: +namespace Step61 { - return 0; -} - -template -class RightHandSide : public Function -{ -public: - RightHandSide() - : Function() - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; -}; + using namespace dealii; -template -double RightHandSide::value(const Point &p, - const unsigned int /*component*/) const -{ - double return_value = 0.0; - return_value = 2 * M_PI * M_PI * sin(M_PI * p[0]) * sin(M_PI * p[1]); - return return_value; -} + // @sect3{The WGDarcyEquation class template} -template -class Solution : public Function -{ -public: - Solution() - : Function(1) - {} - - virtual double value(const Point &p, const unsigned int) const override; -}; + // We will solve for the numerical pressure in the interior and on faces and + // calculate its $L_2$ error of pressure. In the post-processing step, we will + // calculate $L_2$-errors of velocity and flux. + template + class WGDarcyEquation + { + public: + WGDarcyEquation(); + void run(); + + private: + void make_grid(); + void setup_system(); + void assemble_system(); + void solve(); + void postprocess(); + void process_solution(); + void output_results() const; + + Triangulation triangulation; + + AffineConstraints constraints; + + FE_RaviartThomas fe_rt; + DoFHandler dof_handler_rt; + + // The finite element system is used for interior and face solutions. + FESystem fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + // @sect3{Right hand side, boundary values, and exact solution} + + // Next, we define the coefficient matrix $\mathbf{K}$, + // Dirichlet boundary conditions, the right-hand side $f = 2\pi^2 \sin(\pi x) + // \sin(\pi y)$, and the reference solutions $p = \sin(\pi x) \sin(\pi y) $. + // + // The coefficient matrix $\mathbf{K}$ is the identity matrix as a test + // example. + template + class Coefficient : public TensorFunction<2, dim> + { + public: + Coefficient() + : TensorFunction<2, dim>() + {} + + virtual void value_list(const std::vector> &points, + std::vector> &values) const override; + }; + + template + void Coefficient::value_list(const std::vector> &points, + std::vector> & values) const + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); + for (unsigned int p = 0; p < points.size(); ++p) + { + values[p].clear(); + for (unsigned int d = 0; d < dim; ++d) + values[p][d][d] = 1; + } + } -template -double Solution::value(const Point &p, const unsigned int) const -{ - double return_value = 0; - return_value = sin(M_PI * p[0]) * sin(M_PI * p[1]); - return return_value; -} + template + class BoundaryValues : public Function + { + public: + BoundaryValues() + : Function(2) + {} + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + template + double BoundaryValues::value(const Point & /*p*/, + const unsigned int /*component*/) const + { + return 0; + } -template -class Velocity : public TensorFunction<1, dim> -{ -public: - Velocity() - : TensorFunction<1, dim>() - {} + template + class RightHandSide : public Function + { + public: + 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 return_value = 0.0; + return_value = 2 * M_PI * M_PI * sin(M_PI * p[0]) * sin(M_PI * p[1]); + return return_value; + } - virtual Tensor<1, dim> value(const Point &p) const override; -}; + template + class Solution : public Function + { + public: + Solution() + : Function(1) + {} -template -Tensor<1, dim> Velocity::value(const Point &p) const -{ - Tensor<1, dim> return_value; - return_value[0] = -M_PI * cos(M_PI * p[0]) * sin(M_PI * p[1]); - return_value[1] = -M_PI * sin(M_PI * p[0]) * cos(M_PI * p[1]); - return return_value; -} + virtual double value(const Point &p, + const unsigned int) const override; + }; -// @sect3{WGDarcyEquation class implementation} + template + double Solution::value(const Point &p, const unsigned int) const + { + double return_value = 0; + return_value = sin(M_PI * p[0]) * sin(M_PI * p[1]); + return return_value; + } -// @sect4{WGDarcyEquation::WGDarcyEquation} + template + class Velocity : public TensorFunction<1, dim> + { + public: + Velocity() + : TensorFunction<1, dim>() + {} -// In this constructor, we create a finite element space for vector valued -// functions, FE_RaviartThomas. We will need shape functions in -// this space to approximate discrete weak gradients. + virtual Tensor<1, dim> value(const Point &p) const override; + }; -// FESystem defines finite element spaces in the interior and on -// edges of elements. Each of them gets an individual component. Others are the -// same as previous tutorial programs. -template -WGDarcyEquation::WGDarcyEquation() - : fe_rt(0) - , dof_handler_rt(triangulation) - , + template + Tensor<1, dim> Velocity::value(const Point &p) const + { + Tensor<1, dim> return_value; + return_value[0] = -M_PI * cos(M_PI * p[0]) * sin(M_PI * p[1]); + return_value[1] = -M_PI * sin(M_PI * p[0]) * cos(M_PI * p[1]); + return return_value; + } - fe(FE_DGQ(0), 1, FE_FaceQ(0), 1) - , dof_handler(triangulation) + // @sect3{WGDarcyEquation class implementation} -{} + // @sect4{WGDarcyEquation::WGDarcyEquation} -// @sect4{WGDarcyEquation::make_grid} + // In this constructor, we create a finite element space for vector valued + // functions, FE_RaviartThomas. We will need shape functions in + // this space to approximate discrete weak gradients. -// We generate a mesh on the unit square domain and refine it. + // FESystem defines finite element spaces in the interior and on + // edges of elements. Each of them gets an individual component. Others are + // the same as previous tutorial programs. + template + WGDarcyEquation::WGDarcyEquation() + : fe_rt(0) + , dof_handler_rt(triangulation) + , -template -void WGDarcyEquation::make_grid() -{ - GridGenerator::hyper_cube(triangulation, 0, 1); - triangulation.refine_global(1); + fe(FE_DGQ(0), 1, FE_FaceQ(0), 1) + , dof_handler(triangulation) - std::cout << " Number of active cells: " << triangulation.n_active_cells() - << std::endl - << " Total number of cells: " << triangulation.n_cells() - << std::endl; -} + {} -// @sect4{WGDarcyEquation::setup_system} + // @sect4{WGDarcyEquation::make_grid} -// After we create the mesh, we distribute degrees of freedom for the two -// DoFHandler objects. + // We generate a mesh on the unit square domain and refine it. -template -void WGDarcyEquation::setup_system() -{ - dof_handler_rt.distribute_dofs(fe_rt); - dof_handler.distribute_dofs(fe); + template + void WGDarcyEquation::make_grid() + { + GridGenerator::hyper_cube(triangulation, 0, 1); + triangulation.refine_global(1); - std::cout << " Number of flux degrees of freedom: " - << dof_handler_rt.n_dofs() << std::endl; + std::cout << " Number of active cells: " << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " << triangulation.n_cells() + << std::endl; + } - std::cout << " Number of pressure degrees of freedom: " - << dof_handler.n_dofs() << std::endl; + // @sect4{WGDarcyEquation::setup_system} - solution.reinit(dof_handler.n_dofs()); - system_rhs.reinit(dof_handler.n_dofs()); + // After we create the mesh, we distribute degrees of freedom for the two + // DoFHandler objects. + template + void WGDarcyEquation::setup_system() { - constraints.clear(); - FEValuesExtractors::Scalar face(1); - ComponentMask face_pressure_mask = fe.component_mask(face); - VectorTools::interpolate_boundary_values( - dof_handler, 0, BoundaryValues(), constraints, face_pressure_mask); - constraints.close(); - } - - - // In the bilinear form, there is no integration term over faces - // between two neighboring cells, so we can just use - // DoFTools::make_sparsity_pattern to calculate the sparse - // matrix. - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); - sparsity_pattern.copy_from(dsp); + dof_handler_rt.distribute_dofs(fe_rt); + dof_handler.distribute_dofs(fe); - system_matrix.reinit(sparsity_pattern); + std::cout << " Number of flux degrees of freedom: " + << dof_handler_rt.n_dofs() << std::endl; - // solution.reinit(dof_handler.n_dofs()); - // system_rhs.reinit(dof_handler.n_dofs()); -} + std::cout << " Number of pressure degrees of freedom: " + << dof_handler.n_dofs() << std::endl; -// @sect4{WGDarcyEquation::assemble_system} + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); -// First, we create quadrature points and FEValue objects for cells -// and faces. Then we allocate space for all cell matrices and the right-hand -// side vector. The following definitions have been explained in previous -// tutorials. -template -void WGDarcyEquation::assemble_system() -{ - QGauss quadrature_formula(fe_rt.degree + 1); - QGauss face_quadrature_formula(fe_rt.degree + 1); - const RightHandSide right_hand_side; - - // We define objects to evaluate values and - // gradients of shape functions at the quadrature points. - // Since we need shape functions and normal vectors on faces, we need - // FEFaceValues. - FEValues fe_values_rt(fe_rt, - quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - FEValues fe_values(fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face_values(fe, - face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face_values_rt(fe_rt, - face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - - const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - const unsigned int n_q_points = fe_values.get_quadrature().size(); - const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); - const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); - - std::vector local_dof_indices(dofs_per_cell); - - // We will construct these cell matrices to solve for the pressure. - FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); - FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); - FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); - FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_rhs(dofs_per_cell); - Vector cell_solution(dofs_per_cell); - - const Coefficient coefficient; - std::vector> coefficient_values(n_q_points_rt); - - // We need FEValuesExtractors to access the @p interior and - // @p face component of the FESystem shape functions. - const FEValuesExtractors::Vector velocities(0); - const FEValuesExtractors::Scalar interior(0); - const FEValuesExtractors::Scalar face(1); - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - typename DoFHandler::active_cell_iterator cell_rt = - dof_handler_rt.begin_active(); - - // Here, we will calculate cell matrices used to construct the local matrix on - // each cell. We need shape functions for the Raviart-Thomas space as well, so - // we also loop over the corresponding velocity cell iterators. - for (; cell != endc; ++cell, ++cell_rt) { - // On each cell, cell matrices are different, so in every loop, they need - // to be re-computed. - fe_values_rt.reinit(cell_rt); - fe_values.reinit(cell); - coefficient.value_list(fe_values_rt.get_quadrature_points(), - coefficient_values); - - // This cell matrix is the mass matrix for the Raviart-Thomas space. - // Hence, we need to loop over all the quadrature points - // for the velocity FEValues object. - cell_matrix_rt = 0; - for (unsigned int q = 0; q < n_q_points_rt; ++q) - { - for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) - { - const Tensor<1, dim> phi_i_u = - fe_values_rt[velocities].value(i, q); - for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) - { - const Tensor<1, dim> phi_j_u = - fe_values_rt[velocities].value(j, q); - cell_matrix_rt(i, j) += - (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); - } - } - } - // Next we take the inverse of this matrix by using - // gauss_jordan(). It will be used to calculate the - // coefficient matrix later. - cell_matrix_rt.gauss_jordan(); - - // From the introduction, we know that the right hand side - // is the difference between a face integral and a cell integral. - // Here, we approximate the negative of the contribution in the interior. - // Each component of this matrix is the integral of a product between a - // basis function of the polynomial space and the divergence of a basis - // function of the Raviart-Thomas space. These basis functions are defined - // in the interior. - cell_matrix_F = 0; - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const double phi_k_u_div = - fe_values_rt[velocities].divergence(k, q); - cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * - phi_k_u_div * fe_values.JxW(q)); - } - } - } - - // Now, we approximate the integral on faces. - // Each component is the integral of a product between a basis function of - // the polynomial space and the dot product of a basis function of the - // Raviart-Thomas space and the normal vector. So we loop over all the - // faces of the element and obtain the normal vector. - for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; - ++face_n) - { - fe_face_values.reinit(cell, face_n); - fe_face_values_rt.reinit(cell_rt, face_n); - for (unsigned int q = 0; q < n_face_q_points; ++q) - { - const Tensor<1, dim> normal = fe_face_values.normal_vector(q); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const Tensor<1, dim> phi_k_u = - fe_face_values_rt[velocities].value(k, q); - cell_matrix_F(i, k) += - (fe_face_values[face].value(i, q) * (phi_k_u * normal) * - fe_face_values.JxW(q)); - } - } - } - } - - // @p cell_matrix_C is matrix product between the inverse of mass matrix @p cell_matrix_rt and @p cell_matrix_F. - cell_matrix_C = 0; - cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); - - // Element $a_{ij}$ of the local cell matrix $A$ is given by - // $\int_{E} \sum_{k,l} c_{ik} c_{jl} (\mathbf{K} \mathbf{w}_k) \cdot - // \mathbf{w}_l \mathrm{d}x.$ We have calculated coefficients $c$ in the - // previous step. - local_matrix = 0; - for (unsigned int q = 0; q < n_q_points_rt; ++q) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - { - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const Tensor<1, dim> phi_k_u = - fe_values_rt[velocities].value(k, q); - for (unsigned int l = 0; l < dofs_per_cell_rt; ++l) - { - const Tensor<1, dim> phi_l_u = - fe_values_rt[velocities].value(l, q); - local_matrix(i, j) += coefficient_values[q] * - cell_matrix_C[i][k] * phi_k_u * - cell_matrix_C[j][l] * phi_l_u * - fe_values_rt.JxW(q); - } - } - } - } - } + constraints.clear(); + FEValuesExtractors::Scalar face(1); + ComponentMask face_pressure_mask = fe.component_mask(face); + VectorTools::interpolate_boundary_values( + dof_handler, 0, BoundaryValues(), constraints, face_pressure_mask); + constraints.close(); + } - // Next, we calculate the right hand side, $\int_{E} f q \mathrm{d}x$. - cell_rhs = 0; - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - cell_rhs(i) += - (fe_values[interior].value(i, q) * - right_hand_side.value(fe_values.quadrature_point(q)) * - fe_values.JxW(q)); - } - } - // In this part, we distribute components of this local matrix into the - // system matrix and transfer components of the cell right-hand side into - // the system right hand side. - cell->get_dof_indices(local_dof_indices); - constraints.distribute_local_to_global( - local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); - } -} + // In the bilinear form, there is no integration term over faces + // between two neighboring cells, so we can just use + // DoFTools::make_sparsity_pattern to calculate the sparse + // matrix. + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + sparsity_pattern.copy_from(dsp); -// @sect4{WGDarcyEquation::solve} + system_matrix.reinit(sparsity_pattern); -// Solving the system of the Darcy equation. Now, we have pressures in the -// interior and on the faces of all the cells. -template -void WGDarcyEquation::solve() -{ - SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); - SolverCG<> solver(solver_control); - solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); - constraints.distribute(solution); -} + // solution.reinit(dof_handler.n_dofs()); + // system_rhs.reinit(dof_handler.n_dofs()); + } -// @sect4{WGDarcyEquation::process_solution} + // @sect4{WGDarcyEquation::assemble_system} -// This part is to calculate the $L_2$ error of the pressure. -template -void WGDarcyEquation::process_solution() -{ - // Since we have two different spaces for finite elements in interior and on - // faces, if we want to calculate $L_2$ errors in interior, we need degrees of - // freedom only defined in cells. In FESystem, we have two - // components, the first one is for interior, the second one is for skeletons. - // fe.base_element(0) shows we only need degrees of freedom - // defined in cells. - DoFHandler interior_dof_handler(triangulation); - interior_dof_handler.distribute_dofs(fe.base_element(0)); - // We define a vector to extract pressures in cells. - // The size of the vector is the collective number of all degrees of freedom - // in the interior of all the elements. - Vector interior_solution(interior_dof_handler.n_dofs()); + // First, we create quadrature points and FEValue objects for + // cells and faces. Then we allocate space for all cell matrices and the + // right-hand side vector. The following definitions have been explained in + // previous tutorials. + template + void WGDarcyEquation::assemble_system() { - // types::global_dof_index is used to know the global indices - // of degrees of freedom. So here, we get the global indices of local - // degrees of freedom and the global indices of interior degrees of freedom. - std::vector local_dof_indices(fe.dofs_per_cell); - std::vector interior_local_dof_indices( - fe.base_element(0).dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(), - interior_cell = interior_dof_handler.begin_active(); - - // In the loop of all cells and interior of the cell, - // we extract interior solutions from the global solution. - for (; cell != endc; ++cell, ++interior_cell) + QGauss quadrature_formula(fe_rt.degree + 1); + QGauss face_quadrature_formula(fe_rt.degree + 1); + const RightHandSide right_hand_side; + + // We define objects to evaluate values and + // gradients of shape functions at the quadrature points. + // Since we need shape functions and normal vectors on faces, we need + // FEFaceValues. + FEValues fe_values_rt(fe_rt, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values_rt(fe_rt, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + + const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points = fe_values.get_quadrature().size(); + const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); + const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); + + std::vector local_dof_indices(dofs_per_cell); + + // We will construct these cell matrices to solve for the pressure. + FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); + FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + Vector cell_solution(dofs_per_cell); + + const Coefficient coefficient; + std::vector> coefficient_values(n_q_points_rt); + + // We need FEValuesExtractors to access the @p interior and + // @p face component of the FESystem shape functions. + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar interior(0); + const FEValuesExtractors::Scalar face(1); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator cell_rt = + dof_handler_rt.begin_active(); + + // Here, we will calculate cell matrices used to construct the local matrix + // on each cell. We need shape functions for the Raviart-Thomas space as + // well, so we also loop over the corresponding velocity cell iterators. + for (; cell != endc; ++cell, ++cell_rt) { + // On each cell, cell matrices are different, so in every loop, they + // need to be re-computed. + fe_values_rt.reinit(cell_rt); + fe_values.reinit(cell); + coefficient.value_list(fe_values_rt.get_quadrature_points(), + coefficient_values); + + // This cell matrix is the mass matrix for the Raviart-Thomas space. + // Hence, we need to loop over all the quadrature points + // for the velocity FEValues object. + cell_matrix_rt = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) + { + const Tensor<1, dim> phi_i_u = + fe_values_rt[velocities].value(i, q); + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + const Tensor<1, dim> phi_j_u = + fe_values_rt[velocities].value(j, q); + cell_matrix_rt(i, j) += + (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); + } + } + } + // Next we take the inverse of this matrix by using + // gauss_jordan(). It will be used to calculate the + // coefficient matrix later. + cell_matrix_rt.gauss_jordan(); + + // From the introduction, we know that the right hand side + // is the difference between a face integral and a cell integral. + // Here, we approximate the negative of the contribution in the + // interior. Each component of this matrix is the integral of a product + // between a basis function of the polynomial space and the divergence + // of a basis function of the Raviart-Thomas space. These basis + // functions are defined in the interior. + cell_matrix_F = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const double phi_k_u_div = + fe_values_rt[velocities].divergence(k, q); + cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * + phi_k_u_div * fe_values.JxW(q)); + } + } + } + + // Now, we approximate the integral on faces. + // Each component is the integral of a product between a basis function + // of the polynomial space and the dot product of a basis function of + // the Raviart-Thomas space and the normal vector. So we loop over all + // the faces of the element and obtain the normal vector. + for (unsigned int face_n = 0; + face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + cell_matrix_F(i, k) += + (fe_face_values[face].value(i, q) * + (phi_k_u * normal) * fe_face_values.JxW(q)); + } + } + } + } + + // @p cell_matrix_C is matrix product between the inverse of mass matrix @p cell_matrix_rt and @p cell_matrix_F. + cell_matrix_C = 0; + cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); + + // Element $a_{ij}$ of the local cell matrix $A$ is given by + // $\int_{E} \sum_{k,l} c_{ik} c_{jl} (\mathbf{K} \mathbf{w}_k) \cdot + // \mathbf{w}_l \mathrm{d}x.$ We have calculated coefficients $c$ in the + // previous step. + local_matrix = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_values_rt[velocities].value(k, q); + for (unsigned int l = 0; l < dofs_per_cell_rt; ++l) + { + const Tensor<1, dim> phi_l_u = + fe_values_rt[velocities].value(l, q); + local_matrix(i, j) += + coefficient_values[q] * cell_matrix_C[i][k] * + phi_k_u * cell_matrix_C[j][l] * phi_l_u * + fe_values_rt.JxW(q); + } + } + } + } + } + + // Next, we calculate the right hand side, $\int_{E} f q \mathrm{d}x$. + cell_rhs = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + cell_rhs(i) += + (fe_values[interior].value(i, q) * + right_hand_side.value(fe_values.quadrature_point(q)) * + fe_values.JxW(q)); + } + } + + // In this part, we distribute components of this local matrix into the + // system matrix and transfer components of the cell right-hand side + // into the system right hand side. cell->get_dof_indices(local_dof_indices); - interior_cell->get_dof_indices(interior_local_dof_indices); - - for (unsigned int i = 0; i < fe.base_element(0).dofs_per_cell; ++i) - interior_solution(interior_local_dof_indices[i]) = - solution(local_dof_indices[fe.component_to_system_index(0, i)]); + constraints.distribute_local_to_global( + local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } } - // We define a vector that holds the norm of the error on each cell. - // Next, we use VectorTool::integrate_difference - // to compute the error in the $L_2$ norm on each cell. - // Finally, we get the global $L_2$ norm. - Vector difference_per_cell(triangulation.n_active_cells()); - VectorTools::integrate_difference(interior_dof_handler, - interior_solution, - Solution(), - difference_per_cell, - QGauss(fe.degree + 2), - VectorTools::L2_norm); - - const double L2_error = difference_per_cell.l2_norm(); - std::cout << "L2_error_pressure " << L2_error << std::endl; -} - -// @sect4{WGDarcyEquation::postprocess} + // @sect4{WGDarcyEquation::solve} -// After we calculated the numerical pressure, we evaluate $L_2$ errors for the -// velocity on each cell and $L_2$ errors for the flux on faces. + // Solving the system of the Darcy equation. Now, we have pressures in the + // interior and on the faces of all the cells. + template + void WGDarcyEquation::solve() + { + SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); + SolverCG<> solver(solver_control); + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + constraints.distribute(solution); + } -// We are going to evaluate velocities on each cell and calculate the difference -// between numerical and exact velocities. To calculate velocities, we need -// interior and face pressure values of each element, and some other cell -// matrices. + // @sect4{WGDarcyEquation::process_solution} -template -void WGDarcyEquation::postprocess() -{ - QGauss quadrature_formula(fe_rt.degree + 1); - QGauss face_quadrature_formula(fe_rt.degree + 1); - - FEValues fe_values_rt(fe_rt, - quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - FEValues fe_values(fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face_values(fe, - face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - FEFaceValues fe_face_values_rt(fe_rt, - face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); - const unsigned int n_q_points = fe_values.get_quadrature().size(); - const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); - const unsigned int n_face_q_points_rt = - fe_face_values_rt.get_quadrature().size(); - - - std::vector local_dof_indices(dofs_per_cell); - FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); - FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); - FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); - FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); - FullMatrix cell_matrix_D(dofs_per_cell_rt, dofs_per_cell_rt); - FullMatrix cell_matrix_E(dofs_per_cell_rt, dofs_per_cell_rt); - Vector cell_rhs(dofs_per_cell); - Vector cell_solution(dofs_per_cell); - Tensor<1, dim> velocity_cell; - Tensor<1, dim> velocity_face; - Tensor<1, dim> exact_velocity_face; - double L2_err_velocity_cell_sqr_global; - L2_err_velocity_cell_sqr_global = 0; - double L2_err_flux_sqr; - L2_err_flux_sqr = 0; - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - - typename DoFHandler::active_cell_iterator cell_rt = - dof_handler_rt.begin_active(); - - const Coefficient coefficient; - std::vector> coefficient_values(n_q_points_rt); - const FEValuesExtractors::Vector velocities(0); - const FEValuesExtractors::Scalar pressure(dim); - const FEValuesExtractors::Scalar interior(0); - const FEValuesExtractors::Scalar face(1); - - Velocity exact_velocity; - - // In the loop over all cells, we will calculate $L_2$ errors of velocity and - // flux. - - // First, we calculate the $L_2$ velocity error. - // In the introduction, we explained how to calculate the numerical velocity - // on the cell. We need the pressure solution values on each cell, - // coefficients of the Gram matrix and coefficients of the $L_2$ projection. - // We have already calculated the global solution, so we will extract the cell - // solution from the global solution. The coefficients of the Gram matrix have - // been calculated when we assembled the system matrix for the pressures. We - // will do the same way here. For the coefficients of the projection, we do - // matrix multiplication, i.e., the inverse of the Gram matrix times the - // matrix with $(\mathbf{K} \mathbf{w}, \mathbf{w})$ as components. Then, we - // multiply all these coefficients and call them beta. The numerical velocity - // is the product of beta and the basis functions of the Raviart-Thomas space. - for (; cell != endc; ++cell, ++cell_rt) + // This part is to calculate the $L_2$ error of the pressure. + template + void WGDarcyEquation::process_solution() + { + // Since we have two different spaces for finite elements in interior and on + // faces, if we want to calculate $L_2$ errors in interior, we need degrees + // of freedom only defined in cells. In FESystem, we have two + // components, the first one is for interior, the second one is for + // skeletons. fe.base_element(0) shows we only need degrees of + // freedom defined in cells. + DoFHandler interior_dof_handler(triangulation); + interior_dof_handler.distribute_dofs(fe.base_element(0)); + // We define a vector to extract pressures in cells. + // The size of the vector is the collective number of all degrees of freedom + // in the interior of all the elements. + Vector interior_solution(interior_dof_handler.n_dofs()); { - fe_values_rt.reinit(cell_rt); - fe_values.reinit(cell); - coefficient.value_list(fe_values_rt.get_quadrature_points(), - coefficient_values); - - // The component of this cell_matrix_E is the integral of - // $(\mathbf{K} \mathbf{w}, \mathbf{w})$. cell_matrix_rt is - // the Gram matrix. - cell_matrix_E = 0; - cell_matrix_rt = 0; - for (unsigned int q = 0; q < n_q_points_rt; ++q) + // types::global_dof_index is used to know the global indices + // of degrees of freedom. So here, we get the global indices of local + // degrees of freedom and the global indices of interior degrees of + // freedom. + std::vector local_dof_indices(fe.dofs_per_cell); + std::vector interior_local_dof_indices( + fe.base_element(0).dofs_per_cell); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + interior_cell = interior_dof_handler.begin_active(); + + // In the loop of all cells and interior of the cell, + // we extract interior solutions from the global solution. + for (; cell != endc; ++cell, ++interior_cell) { - for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) - { - const Tensor<1, dim> phi_i_u = - fe_values_rt[velocities].value(i, q); - - for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) - { - const Tensor<1, dim> phi_j_u = - fe_values_rt[velocities].value(j, q); - - cell_matrix_E(i, j) += (coefficient_values[q] * phi_j_u * - phi_i_u * fe_values_rt.JxW(q)); - cell_matrix_rt(i, j) += - (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); - } - } + cell->get_dof_indices(local_dof_indices); + interior_cell->get_dof_indices(interior_local_dof_indices); + + for (unsigned int i = 0; i < fe.base_element(0).dofs_per_cell; ++i) + interior_solution(interior_local_dof_indices[i]) = + solution(local_dof_indices[fe.component_to_system_index(0, i)]); } + } - // We take the inverse of the Gram matrix, take matrix multiplication and - // get the matrix with coefficients of projection. - cell_matrix_D = 0; - cell_matrix_rt.gauss_jordan(); - cell_matrix_rt.mmult(cell_matrix_D, cell_matrix_E); + // We define a vector that holds the norm of the error on each cell. + // Next, we use VectorTool::integrate_difference + // to compute the error in the $L_2$ norm on each cell. + // Finally, we get the global $L_2$ norm. + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(interior_dof_handler, + interior_solution, + Solution(), + difference_per_cell, + QGauss(fe.degree + 2), + VectorTools::L2_norm); + + const double L2_error = difference_per_cell.l2_norm(); + std::cout << "L2_error_pressure " << L2_error << std::endl; + } - // This cell matrix will be used to calculate the coefficients of the Gram - // matrix. This part is the same as the part in evaluating pressure. - cell_matrix_F = 0; - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const double phi_k_u_div = - fe_values_rt[velocities].divergence(k, q); - cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * - phi_k_u_div * fe_values.JxW(q)); - } - } - } + // @sect4{WGDarcyEquation::postprocess} - for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; - ++face_n) - { - fe_face_values.reinit(cell, face_n); - fe_face_values_rt.reinit(cell_rt, face_n); - for (unsigned int q = 0; q < n_face_q_points; ++q) - { - const Tensor<1, dim> normal = fe_face_values.normal_vector(q); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const Tensor<1, dim> phi_k_u = - fe_face_values_rt[velocities].value(k, q); - cell_matrix_F(i, k) += - (fe_face_values[face].value(i, q) * (phi_k_u * normal) * - fe_face_values.JxW(q)); - } - } - } - } - cell_matrix_C = 0; - cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); + // After we calculated the numerical pressure, we evaluate $L_2$ errors for + // the velocity on each cell and $L_2$ errors for the flux on faces. - // This is to extract pressure values of the element. - cell->get_dof_indices(local_dof_indices); - cell_solution = 0; - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - cell_solution(i) = solution(local_dof_indices[i]); - } + // We are going to evaluate velocities on each cell and calculate the + // difference between numerical and exact velocities. To calculate velocities, + // we need interior and face pressure values of each element, and some other + // cell matrices. - // From previous calculations we obtained all the coefficients needed to - // calculate beta. - Vector beta(dofs_per_cell_rt); - beta = 0; - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - beta(k) += -(cell_solution(i) * cell_matrix_C(i, j) * - cell_matrix_D(k, j)); - } - } - } + template + void WGDarcyEquation::postprocess() + { + QGauss quadrature_formula(fe_rt.degree + 1); + QGauss face_quadrature_formula(fe_rt.degree + 1); + + FEValues fe_values_rt(fe_rt, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values_rt(fe_rt, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); + const unsigned int n_q_points = fe_values.get_quadrature().size(); + const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); + const unsigned int n_face_q_points_rt = + fe_face_values_rt.get_quadrature().size(); + + + std::vector local_dof_indices(dofs_per_cell); + FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); + FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix_D(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_E(dofs_per_cell_rt, dofs_per_cell_rt); + Vector cell_rhs(dofs_per_cell); + Vector cell_solution(dofs_per_cell); + Tensor<1, dim> velocity_cell; + Tensor<1, dim> velocity_face; + Tensor<1, dim> exact_velocity_face; + double L2_err_velocity_cell_sqr_global; + L2_err_velocity_cell_sqr_global = 0; + double L2_err_flux_sqr; + L2_err_flux_sqr = 0; + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + + typename DoFHandler::active_cell_iterator cell_rt = + dof_handler_rt.begin_active(); + + const Coefficient coefficient; + std::vector> coefficient_values(n_q_points_rt); + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + const FEValuesExtractors::Scalar interior(0); + const FEValuesExtractors::Scalar face(1); + + Velocity exact_velocity; + + // In the loop over all cells, we will calculate $L_2$ errors of velocity + // and flux. + + // First, we calculate the $L_2$ velocity error. + // In the introduction, we explained how to calculate the numerical velocity + // on the cell. We need the pressure solution values on each cell, + // coefficients of the Gram matrix and coefficients of the $L_2$ projection. + // We have already calculated the global solution, so we will extract the + // cell solution from the global solution. The coefficients of the Gram + // matrix have been calculated when we assembled the system matrix for the + // pressures. We will do the same way here. For the coefficients of the + // projection, we do matrix multiplication, i.e., the inverse of the Gram + // matrix times the matrix with $(\mathbf{K} \mathbf{w}, \mathbf{w})$ as + // components. Then, we multiply all these coefficients and call them beta. + // The numerical velocity is the product of beta and the basis functions of + // the Raviart-Thomas space. + for (; cell != endc; ++cell, ++cell_rt) + { + fe_values_rt.reinit(cell_rt); + fe_values.reinit(cell); + coefficient.value_list(fe_values_rt.get_quadrature_points(), + coefficient_values); + + // The component of this cell_matrix_E is the integral of + // $(\mathbf{K} \mathbf{w}, \mathbf{w})$. cell_matrix_rt is + // the Gram matrix. + cell_matrix_E = 0; + cell_matrix_rt = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) + { + const Tensor<1, dim> phi_i_u = + fe_values_rt[velocities].value(i, q); + + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + const Tensor<1, dim> phi_j_u = + fe_values_rt[velocities].value(j, q); + + cell_matrix_E(i, j) += (coefficient_values[q] * phi_j_u * + phi_i_u * fe_values_rt.JxW(q)); + cell_matrix_rt(i, j) += + (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); + } + } + } + + // We take the inverse of the Gram matrix, take matrix multiplication + // and get the matrix with coefficients of projection. + cell_matrix_D = 0; + cell_matrix_rt.gauss_jordan(); + cell_matrix_rt.mmult(cell_matrix_D, cell_matrix_E); + + // This cell matrix will be used to calculate the coefficients of the + // Gram matrix. This part is the same as the part in evaluating + // pressure. + cell_matrix_F = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const double phi_k_u_div = + fe_values_rt[velocities].divergence(k, q); + cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * + phi_k_u_div * fe_values.JxW(q)); + } + } + } + + for (unsigned int face_n = 0; + face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + cell_matrix_F(i, k) += + (fe_face_values[face].value(i, q) * + (phi_k_u * normal) * fe_face_values.JxW(q)); + } + } + } + } + cell_matrix_C = 0; + cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); + + // This is to extract pressure values of the element. + cell->get_dof_indices(local_dof_indices); + cell_solution = 0; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + cell_solution(i) = solution(local_dof_indices[i]); + } + + // From previous calculations we obtained all the coefficients needed to + // calculate beta. + Vector beta(dofs_per_cell_rt); + beta = 0; + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + beta(k) += -(cell_solution(i) * cell_matrix_C(i, j) * + cell_matrix_D(k, j)); + } + } + } + + // Now, we can calculate the numerical velocity at each quadrature point + // and compute the $L_2$ error on each cell. + double L2_err_velocity_cell_sqr_local; + double difference_velocity_cell_sqr; + L2_err_velocity_cell_sqr_local = 0; + velocity_cell = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + difference_velocity_cell_sqr = 0; + velocity_cell = 0; + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_values_rt[velocities].value(k, q); + velocity_cell += beta(k) * phi_k_u; + } + difference_velocity_cell_sqr = + (velocity_cell - + exact_velocity.value(fe_values_rt.quadrature_point(q))) * + (velocity_cell - + exact_velocity.value(fe_values_rt.quadrature_point(q))); + L2_err_velocity_cell_sqr_local += + difference_velocity_cell_sqr * fe_values_rt.JxW(q); + } + + L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local; + + // For reconstructing the flux we need the size of cells and faces. + // Since fluxes are calculated on faces, we have the loop over all four + // faces of each cell. To calculate face velocity, we use the + // coefficient beta we have calculated previously. Then, we calculate + // the squared velocity error in normal direction. Finally, we calculate + // $L_2$ flux error on the cell and add it to the global error. + double difference_velocity_face_sqr; + double L2_err_flux_face_sqr_local; + double err_flux_each_face; + double err_flux_face; + L2_err_flux_face_sqr_local = 0; + err_flux_face = 0; + const double cell_area = cell->measure(); + for (unsigned int face_n = 0; + face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + const double face_length = cell->face(face_n)->measure(); + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + L2_err_flux_face_sqr_local = 0; + err_flux_each_face = 0; + for (unsigned int q = 0; q < n_face_q_points_rt; ++q) + { + difference_velocity_face_sqr = 0; + velocity_face = 0; + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + velocity_face += beta(k) * phi_k_u; + } + exact_velocity_face = + exact_velocity.value(fe_face_values_rt.quadrature_point(q)); + difference_velocity_face_sqr = + (velocity_face * normal - exact_velocity_face * normal) * + (velocity_face * normal - exact_velocity_face * normal); + L2_err_flux_face_sqr_local += + difference_velocity_face_sqr * fe_face_values_rt.JxW(q); + } + err_flux_each_face = + L2_err_flux_face_sqr_local / (face_length) * (cell_area); + err_flux_face += err_flux_each_face; + } + L2_err_flux_sqr += err_flux_face; + } - // Now, we can calculate the numerical velocity at each quadrature point - // and compute the $L_2$ error on each cell. - double L2_err_velocity_cell_sqr_local; - double difference_velocity_cell_sqr; - L2_err_velocity_cell_sqr_local = 0; - velocity_cell = 0; - for (unsigned int q = 0; q < n_q_points_rt; ++q) - { - difference_velocity_cell_sqr = 0; - velocity_cell = 0; - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const Tensor<1, dim> phi_k_u = - fe_values_rt[velocities].value(k, q); - velocity_cell += beta(k) * phi_k_u; - } - difference_velocity_cell_sqr = - (velocity_cell - - exact_velocity.value(fe_values_rt.quadrature_point(q))) * - (velocity_cell - - exact_velocity.value(fe_values_rt.quadrature_point(q))); - L2_err_velocity_cell_sqr_local += - difference_velocity_cell_sqr * fe_values_rt.JxW(q); - } + // After adding up errors over all cells, we take square root and get the + // $L_2$ errors of velocity and flux. + const double L2_err_velocity_cell = + std::sqrt(L2_err_velocity_cell_sqr_global); + std::cout << "L2_error_vel " << L2_err_velocity_cell << std::endl; + const double L2_err_flux_face = std::sqrt(L2_err_flux_sqr); + std::cout << "L2_error_flux " << L2_err_flux_face << std::endl; + } - L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local; - - // For reconstructing the flux we need the size of cells and faces. Since - // fluxes are calculated on faces, we have the loop over all four faces of - // each cell. To calculate face velocity, we use the coefficient beta we - // have calculated previously. Then, we calculate the squared velocity - // error in normal direction. Finally, we calculate $L_2$ flux error on - // the cell and add it to the global error. - double difference_velocity_face_sqr; - double L2_err_flux_face_sqr_local; - double err_flux_each_face; - double err_flux_face; - L2_err_flux_face_sqr_local = 0; - err_flux_face = 0; - const double cell_area = cell->measure(); - for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; - ++face_n) - { - const double face_length = cell->face(face_n)->measure(); - fe_face_values.reinit(cell, face_n); - fe_face_values_rt.reinit(cell_rt, face_n); - L2_err_flux_face_sqr_local = 0; - err_flux_each_face = 0; - for (unsigned int q = 0; q < n_face_q_points_rt; ++q) - { - difference_velocity_face_sqr = 0; - velocity_face = 0; - const Tensor<1, dim> normal = fe_face_values.normal_vector(q); - for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) - { - const Tensor<1, dim> phi_k_u = - fe_face_values_rt[velocities].value(k, q); - velocity_face += beta(k) * phi_k_u; - } - exact_velocity_face = - exact_velocity.value(fe_face_values_rt.quadrature_point(q)); - difference_velocity_face_sqr = - (velocity_face * normal - exact_velocity_face * normal) * - (velocity_face * normal - exact_velocity_face * normal); - L2_err_flux_face_sqr_local += - difference_velocity_face_sqr * fe_face_values_rt.JxW(q); - } - err_flux_each_face = - L2_err_flux_face_sqr_local / (face_length) * (cell_area); - err_flux_face += err_flux_each_face; - } - L2_err_flux_sqr += err_flux_face; - } - // After adding up errors over all cells, we take square root and get the - // $L_2$ errors of velocity and flux. - const double L2_err_velocity_cell = - std::sqrt(L2_err_velocity_cell_sqr_global); - std::cout << "L2_error_vel " << L2_err_velocity_cell << std::endl; - const double L2_err_flux_face = std::sqrt(L2_err_flux_sqr); - std::cout << "L2_error_flux " << L2_err_flux_face << std::endl; -} + // @sect4{WGDarcyEquation::output_results} + // We have 2 sets of results to output: the interior solution + // and the skeleton solution. We use DataOut to visualize + // interior results. The graphical output for the skeleton results is done by + // using the DataOutFaces class. + template + void WGDarcyEquation::output_results() const + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "Pressure_Interior"); + data_out.build_patches(fe.degree); + std::ofstream output("Pressure_Interior.vtk"); + data_out.write_vtk(output); + + DataOutFaces data_out_face(false); + std::vector + face_component_type(2, DataComponentInterpretation::component_is_scalar); + data_out_face.add_data_vector(dof_handler, + solution, + "Pressure_Edge", + face_component_type); + data_out_face.build_patches(fe.degree); + std::ofstream face_output("Pressure_Edge.vtk"); + data_out_face.write_vtk(face_output); + } -// @sect4{WGDarcyEquation::output_results} -// We have 2 sets of results to output: the interior solution -// and the skeleton solution. We use DataOut to visualize interior -// results. The graphical output for the skeleton results is done by using the -// DataOutFaces class. -template -void WGDarcyEquation::output_results() const -{ - DataOut data_out; - data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "Pressure_Interior"); - data_out.build_patches(fe.degree); - std::ofstream output("Pressure_Interior.vtk"); - data_out.write_vtk(output); - - DataOutFaces data_out_face(false); - std::vector - face_component_type(2, DataComponentInterpretation::component_is_scalar); - data_out_face.add_data_vector(dof_handler, - solution, - "Pressure_Edge", - face_component_type); - data_out_face.build_patches(fe.degree); - std::ofstream face_output("Pressure_Edge.vtk"); - data_out_face.write_vtk(face_output); -} + // @sect4{WGDarcyEquation::run} + // This is the final function of the main class. It calls the other functions + // of our class. + template + void WGDarcyEquation::run() + { + std::cout << "Solving problem in " << dim << " space dimensions." + << std::endl; + make_grid(); + setup_system(); + assemble_system(); + solve(); + process_solution(); + postprocess(); + output_results(); + } -// @sect4{WGDarcyEquation::run} +} // namespace Step61 -// This is the final function of the main class. It calls the other functions of -// our class. -template -void WGDarcyEquation::run() -{ - std::cout << "Solving problem in " << dim << " space dimensions." - << std::endl; - make_grid(); - setup_system(); - assemble_system(); - solve(); - process_solution(); - postprocess(); - output_results(); -} // @sect3{The main function} @@ -914,9 +931,9 @@ int main() { try { - deallog.depth_console(2); - WGDarcyEquation<2> WGDarcyEquationTest; - WGDarcyEquationTest.run(); + dealii::deallog.depth_console(2); + Step61::WGDarcyEquation<2> wg_darcy; + wg_darcy.run(); } catch (std::exception &exc) { @@ -929,6 +946,7 @@ int main() << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; + return 1; } catch (...) { @@ -940,7 +958,7 @@ int main() << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - throw; + return 1; } return 0; -- 2.39.5