From: Jiaqi Zhang Date: Mon, 31 May 2021 03:09:15 +0000 (-0400) Subject: feinterface: stokes test X-Git-Tag: v9.4.0-rc1~1235^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12370%2Fhead;p=dealii.git feinterface: stokes test --- diff --git a/tests/feinterface/stokes.cc b/tests/feinterface/stokes.cc new file mode 100644 index 0000000000..8c6783e622 --- /dev/null +++ b/tests/feinterface/stokes.cc @@ -0,0 +1,921 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2020 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * Author: + * Timo Heister, Clemson University + */ + +// test FEInterfaceValues for a DG Stokes problem. + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + +namespace StokesTests +{ + using namespace dealii; + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + std::array values; + std::array cell_indices; + }; + + + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + double value; + unsigned int cell_index; + + + + template + void + reinit(const Iterator &cell, const unsigned int dofs_per_cell) + { + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } + }; + + + + template + inline void + copy(const CopyData & c, + const AffineConstraints &constraints, + MatrixType & system_matrix, + VectorType & system_rhs) + { + constraints.distribute_local_to_global(c.cell_matrix, + c.cell_rhs, + c.local_dof_indices, + system_matrix, + system_rhs); + for (auto &cdf : c.face_data) + { + const unsigned int dofs_per_cell = cdf.joint_dof_indices.size(); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int k = 0; k < dofs_per_cell; ++k) + system_matrix.add(cdf.joint_dof_indices[i], + cdf.joint_dof_indices[k], + cdf.cell_matrix(i, k)); + } + } + + // @sect3{Functions for Solution and Righthand side} + // + // The class Solution is used to define the boundary conditions and to + // compute errors of the numerical solution. Note that we need to define the + // values and gradients in order to compute L2 and H1 errors. Here we + // decided to separate the implementations for 2d and 3d using template + // specialization. + // + // Note that the first dim components are the velocity components + // and the last is the pressure. + template + class Solution : public Function + { + public: + Solution() + : Function(dim + 1) + {} + virtual double + value(const Point &p, const unsigned int component = 0) const; + virtual Tensor<1, dim> + gradient(const Point &p, const unsigned int component = 0) const; + }; + + template <> + double + Solution<2>::value(const Point<2> &p, const unsigned int component) const + { + using numbers::PI; + const double x = p(0); + const double y = p(1); + // zero on BD's + if (component == 0) + return PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y); + if (component == 1) + return -PI * sin(PI * y) * sin(PI * y) * sin(2.0 * PI * x); + if (component == 2) + return cos(PI * x) * sin(PI * y); + + return 0; + } + + template <> + double + Solution<3>::value(const Point<3> &p, const unsigned int component) const + { + Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + const double z = p(2); + + if (component == 0) + return 2. * PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y) * + sin(2.0 * PI * z); + if (component == 1) + return -PI * sin(PI * y) * sin(PI * y) * sin(2.0 * PI * x) * + sin(2.0 * PI * z); + if (component == 2) + return -PI * sin(PI * z) * sin(PI * z) * sin(2.0 * PI * x) * + sin(2.0 * PI * y); + if (component == 3) + return sin(PI * x) * cos(PI * y) * sin(PI * z); + + return 0; + } + + // Note that for the gradient we need to return a Tensor<1,dim> + template <> + Tensor<1, 2> + Solution<2>::gradient(const Point<2> &p, const unsigned int component) const + { + Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + + Tensor<1, 2> return_value; + if (component == 0) + { + return_value[0] = PI * PI * sin(2.0 * PI * y) * sin(2.0 * PI * x); + return_value[1] = + 2.0 * PI * PI * sin(PI * x) * sin(PI * x) * cos(2.0 * PI * y); + } + else if (component == 1) + { + return_value[0] = + -2.0 * PI * PI * sin(PI * y) * sin(PI * y) * cos(2.0 * PI * x); + return_value[1] = -PI * PI * sin(2.0 * PI * y) * sin(2.0 * PI * x); + } + else if (component == 2) + { + return_value[0] = PI * cos(PI * x) * cos(PI * y); + return_value[1] = -PI * sin(PI * x) * sin(PI * y); + } + + return return_value; + } + + template <> + Tensor<1, 3> + Solution<3>::gradient(const Point<3> &p, const unsigned int component) const + { + Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + const double z = p(2); + + Tensor<1, 3> return_value; + if (component == 0) + { + return_value[0] = + 2 * PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z); + return_value[1] = 4 * PI * PI * sin(PI * x) * sin(PI * x) * + cos(2 * PI * y) * sin(2 * PI * z); + return_value[2] = 4 * PI * PI * sin(PI * x) * sin(PI * x) * + cos(2 * PI * z) * sin(2 * PI * y); + } + else if (component == 1) + { + return_value[0] = -2 * PI * PI * sin(PI * y) * sin(PI * y) * + cos(2 * PI * x) * sin(2 * PI * z); + return_value[1] = + -PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z); + return_value[2] = -2 * PI * PI * sin(PI * y) * sin(PI * y) * + cos(2 * PI * z) * sin(2 * PI * x); + } + else if (component == 2) + { + return_value[0] = -2 * PI * PI * sin(PI * z) * sin(PI * z) * + cos(2 * PI * x) * sin(2 * PI * y); + return_value[1] = -2 * PI * PI * sin(PI * z) * sin(PI * z) * + cos(2 * PI * y) * sin(2 * PI * x); + return_value[2] = + -PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z); + } + else if (component == 3) + { + return_value[0] = PI * cos(PI * x) * cos(PI * y) * sin(PI * z); + return_value[1] = -PI * sin(PI * x) * sin(PI * y) * sin(PI * z); + return_value[2] = PI * sin(PI * x) * cos(PI * y) * cos(PI * z); + } + + return return_value; + } + + + + // Implementation of $f$. See the introduction for more information. + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function(dim + 1) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; + }; + + template <> + double + RightHandSide<2>::value(const Point<2> &p, const unsigned int component) const + { + Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); + + using numbers::PI; + double x = p(0); + double y = p(1); + double nu = 1.0; + + // RHS for 0 BD's + if (component == 0) + return -nu * 2.0 * PI * PI * PI * + (-2.0 * sin(PI * x) * sin(PI * x) + cos(2. * PI * x)) * + sin(2.0 * PI * y) - + PI * sin(PI * x) * sin(PI * y); + if (component == 1) + return nu * 2.0 * PI * PI * PI * (2.0 * cos(2.0 * PI * y) - 1) * + sin(2.0 * PI * x) + + PI * cos(PI * x) * cos(PI * y); + if (component == 2) + return 0.0; + + return 0.0; + } + + template <> + double + RightHandSide<3>::value(const Point<3> &p, const unsigned int component) const + { + Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + double x = p(0); + double y = p(1); + double z = p(2); + + if (component == 0) + return 4. * PI * PI * PI * + (4. * sin(PI * x) * sin(PI * x) - cos(2. * PI * x)) * + sin(2. * PI * y) * sin(2. * PI * z) + + PI * cos(PI * x) * cos(PI * y) * sin(PI * z); + if (component == 1) + return -2. * PI * PI * PI * + (4. * sin(PI * y) * sin(PI * y) - cos(2. * PI * y)) * + sin(2. * PI * x) * sin(2. * PI * z) + + PI * (-1) * sin(PI * y) * sin(PI * x) * sin(PI * z); + if (component == 2) + return -2. * PI * PI * PI * + (4. * sin(PI * z) * sin(PI * z) - cos(2. * PI * z)) * + sin(2. * PI * x) * sin(2. * PI * y) + + PI * cos(PI * z) * sin(PI * x) * cos(PI * y); + if (component == 3) + return 0.0; + + return 0.0; + } + + + + // @sect3{The StokesProblem class} + // + // This is the main class of the problem. + template + class StokesProblem + { + public: + StokesProblem(FiniteElement &fe, const unsigned int pressure_degree); + void + run(); + + private: + void + setup_dofs(); + void + assemble_system_mesh_loop(); + void + solve(); + void + compute_errors(unsigned int k); + + const unsigned int pressure_degree; + + Triangulation triangulation; + FiniteElement &fe; + DoFHandler dof_handler; + + AffineConstraints<> constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + SparseMatrix pressure_mass_matrix; + + BlockVector solution; + BlockVector system_rhs; + + double last_l2_error; + double last_H1_error; + double last_Hdiv_error1; + double last_Hdiv_error2; + }; + + + + template + StokesProblem::StokesProblem(FiniteElement &fe, + const unsigned int pressure_degree) + : pressure_degree(pressure_degree) + , triangulation(Triangulation::maximum_smoothing) + , fe(fe) + , dof_handler(triangulation) + {} + + // @sect4{StokesProblem::setup_dofs} + + // This function sets up the DoFHandler, matrices, vectors, and Multigrid + // structures (if needed). + template + void + StokesProblem::setup_dofs() + { + system_matrix.clear(); + pressure_mass_matrix.clear(); + + // The main DoFHandler only needs active DoFs, so we are not calling + // distribute_mg_dofs() here + dof_handler.distribute_dofs(fe); + + // This block structure separates the dim velocity components from + // the pressure component (used for reordering). Note that we have + // 2 instead of dim+1 blocks like in step-22, because our FESystem + // is nested and the dim velocity components appear as one block. + std::vector block_component(2); + block_component[0] = 0; + block_component[1] = 1; + + // Velocities start at component 0: + const FEValuesExtractors::Vector velocities(0); + + // ILU behaves better if we apply a reordering to reduce fillin. There + // is no advantage in doing this for the other solvers. + DoFRenumbering::Cuthill_McKee(dof_handler); + + + // This ensures that all velocities DoFs are enumerated before the + // pressure unknowns. This allows us to use blocks for vectors and + // matrices and allows us to get the same DoF numbering for + // dof_handler and velocity_dof_handler. + DoFRenumbering::block_wise(dof_handler); + + std::vector dofs_per_block(2); + DoFTools::count_dofs_per_block(dof_handler, + dofs_per_block, + block_component); + const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; + + deallog << "\tNumber of active cells: " << triangulation.n_active_cells() + << std::endl + << "\tNumber of degrees of freedom: " << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p << ')' << std::endl; + + { + constraints.reinit(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Solution(), + constraints, + fe.component_mask(velocities)); + constraints.close(); + } + + { + BlockDynamicSparsityPattern csp(dofs_per_block, dofs_per_block); + DoFTools::make_flux_sparsity_pattern(dof_handler, csp, constraints); + sparsity_pattern.copy_from(csp); + } + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dofs_per_block); + system_rhs.reinit(dofs_per_block); + } + + + + template + void + StokesProblem::assemble_system_mesh_loop() + { + system_matrix = 0; + system_rhs = 0; + + typedef decltype(dof_handler.begin_active()) Iterator; + const RightHandSide rhs_function; + const Solution boundary_function; + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + const double nu = 1.0; + + auto penalty_parameter = [](const double degree, + const double extent1, + const double extent2) -> double { + return 4.0 * degree * (degree + 1.0) * 0.5 * + (1.0 / extent1 + 1.0 / extent2); + }; + + auto cell_worker = [&](const Iterator & cell, + MeshWorker::ScratchData &scratch_data, + CopyData & copy_data) { + const FEValues &fe_v = scratch_data.reinit(cell); + + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; + const unsigned int n_q_points = fe_v.get_quadrature().size(); + + copy_data.reinit(cell, dofs_per_cell); + + const std::vector &JxW = fe_v.get_JxW_values(); + + const double nu = 1.0; + std::vector> rhs_values(n_q_points, + Vector(dim + 1)); + rhs_function.vector_value_list(fe_v.get_quadrature_points(), rhs_values); + Tensor<1, dim> force_f; + + for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point) + { + for (unsigned int d = 0; d < dim; ++d) + force_f[d] = rhs_values[point](d); + for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + ( + // nu \nabla v : \nabla u + nu * scalar_product(fe_v[velocities].gradient(i, point), + fe_v[velocities].gradient(j, point)) + // -q, div u + - fe_v[pressure].value(i, point) * + fe_v[velocities].divergence(j, point) + // -p, div v + - fe_v[pressure].value(j, point) * + fe_v[velocities].divergence(i, point) + // p,q + + fe_v[pressure].value(j, point) * + fe_v[pressure].value(i, point)) * + JxW[point]; + + copy_data.cell_rhs(i) += + // f,v + (force_f * fe_v[velocities].value(i, point)) * JxW[point]; + } + } + }; + + auto boundary_worker = [&](const Iterator & cell, + const unsigned int & face_no, + MeshWorker::ScratchData &scratch_data, + CopyData & copy_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + + const auto &q_points = fe_fv.get_quadrature_points(); + + const std::vector & JxW = fe_fv.get_JxW_values(); + const std::vector> &normals = fe_fv.get_normal_vectors(); + + std::vector> g_values(q_points.size(), + Vector(dim + 1)); + boundary_function.vector_value_list(q_points, g_values); + Tensor<1, dim> g; + + const double degree = + std::max(1.0, static_cast(fe_fv.get_fe().degree)); + const double extent1 = cell->extent_in_direction( + GeometryInfo::unit_normal_direction[face_no]); + const double penalty = penalty_parameter(degree, extent1, extent1); + + for (unsigned int point = 0; point < q_points.size(); ++point) + { + for (unsigned int d = 0; d < dim; ++d) + g[d] = g_values[point](d); + + for (unsigned int i = 0; i < fe_fv.dofs_per_cell; ++i) + for (unsigned int j = 0; j < fe_fv.dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + ( + // - nu (\nabla u n) . v + -nu * + ((fe_fv[velocities].gradient(j, point) * normals[point]) * + fe_fv[velocities].value(i, point)) + + // - nu u . (\nabla v n) // NIPG: use + + - + nu * (fe_fv[velocities].value(j, point) * + (fe_fv[velocities].gradient(i, point) * normals[point])) + + // + nu * penalty u . v + + nu * penalty * + (fe_fv[velocities].value(j, point) * + fe_fv[velocities].value(i, point)) + + // p (v.n) + + fe_fv[pressure].value(j, point) * + scalar_product(fe_fv[velocities].value(i, point), + normals[point]) + + // q (u.n) + + fe_fv[pressure].value(i, point) * + scalar_product(fe_fv[velocities].value(j, point), + normals[point]) + + ) * + JxW[point]; + + for (unsigned int i = 0; i < fe_fv.dofs_per_cell; ++i) + copy_data.cell_rhs(i) += + ( + // -nu g . (\nabla v n) // NIPG: use + + -nu * scalar_product(g, + (fe_fv[velocities].gradient(i, point) * + normals[point])) + + // +nu penalty g . v + + nu * penalty * + scalar_product(g, fe_fv[velocities].value(i, point)) + + // q (g.n) (weak normal component of boundary condition) + + fe_fv[pressure].value(i, point) * + scalar_product(g, normals[point])) * + JxW[point]; + } + }; + + auto face_worker = [&](const Iterator & cell, + const unsigned int & f, + const unsigned int & sf, + const Iterator & ncell, + const unsigned int & nf, + const unsigned int & nsf, + MeshWorker::ScratchData &scratch_data, + CopyData & copy_data) { + const FEInterfaceValues &fe_fv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + FEInterfaceViews::Scalar interface_scalar(fe_fv, pressure.component); + FEInterfaceViews::Vector interface_vector( + fe_fv, velocities.first_vector_component); + + copy_data.face_data.emplace_back(); + CopyDataFace & copy_data_face = copy_data.face_data.back(); + const unsigned int dofs_per_cell = fe_fv.n_current_interface_dofs(); + + copy_data_face.joint_dof_indices = fe_fv.get_interface_dof_indices(); + copy_data_face.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + + const std::vector & JxW = fe_fv.get_JxW_values(); + const std::vector> &normals = fe_fv.get_normal_vectors(); + const auto &q_points = fe_fv.get_quadrature_points(); + + double nu = 1.0; + const double degree = + std::max(1.0, static_cast(fe_fv.get_fe().degree)); + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = penalty_parameter(degree, extent1, extent2); + + for (unsigned int point = 0; point < q_points.size(); ++point) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data_face.cell_matrix(i, j) += + ( + // - nu {\nabla u}n . [v] (consistency) + -nu * + (fe_fv[velocities].average_gradient(j, point) * + normals[point]) * + fe_fv[velocities].jump(i, point) + + // - nu [u] . {\nabla v}n (symmetry) // NIPG: use + + - nu * fe_fv[velocities].jump(j, point) * + (fe_fv[velocities].average_gradient(i, point) * + normals[point]) + + // nu sigma [u].[v] (penalty) + + nu * penalty * + scalar_product(fe_fv[velocities].jump(j, point), + fe_fv[velocities].jump(i, point)) + + // {p} ([v].n) + + fe_fv[pressure].average(j, point) * + scalar_product(fe_fv[velocities].jump(i, point), + normals[point]) + + // {q} ([u].n) + + fe_fv[pressure].average(i, point) * + scalar_product(fe_fv[velocities].jump(j, point), + normals[point])) * + JxW[point]; + } + }; + + auto copier = [&](const CopyData &c) { + copy(c, constraints, system_matrix, system_rhs); + }; + + const unsigned int n_gauss_points = pressure_degree + 2; + const UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + const UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | + update_normal_vectors | update_JxW_values; + + const QGauss quadrature(n_gauss_points); + const QGauss face_quadrature(n_gauss_points); + + static MappingQ1 mapping; + MeshWorker::ScratchData scratch_data( + mapping, fe, quadrature, cell_flags, face_quadrature, face_flags); + CopyData cd; + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + cd, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_boundary_faces | + MeshWorker::assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + pressure_mass_matrix.reinit(sparsity_pattern.block(1, 1)); + pressure_mass_matrix.copy_from(system_matrix.block(1, 1)); + system_matrix.block(1, 1) = 0; + } + + + + template + void + StokesProblem::solve() + { + solution = 0; + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + + constraints.distribute(solution); + } + + + + // @sect4{StokesProblem::process_solution} + + // This function computes the L2 and H1 errors of the solution. For this, + // we need to make sure the pressure has mean zero. + template + void + StokesProblem::compute_errors(unsigned int k) + { + const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), + dim + 1); + const ComponentSelectFunction pressure_mask(dim, 1.0, dim + 1); + + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 3), + VectorTools::L2_norm, + &velocity_mask); + + const double Velocity_L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 3), + VectorTools::H1_norm, + &velocity_mask); + + const double Velocity_H1_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference(dof_handler, + solution, + Functions::ZeroFunction(dim + 1), + difference_per_cell, + QGauss(pressure_degree + 3), + VectorTools::Hdiv_seminorm, + &velocity_mask); + + const double Velocity_Hdiv_error1 = difference_per_cell.l2_norm(); + + static double last_Pressure_L2_error = 0; + + const double mean_pressure = VectorTools::compute_mean_value( + dof_handler, QGauss(pressure_degree + 3), solution, dim); + solution.block(1).add(-mean_pressure); + + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 3), + VectorTools::L2_norm, + &pressure_mask); + const double Pressure_L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + + deallog << " At " << k + 1 << "th mesh" << std::endl + << " L2 error: " << std::setw(12) << Velocity_L2_error + << std::setw(0) << " L2_Conv_rate: " << std::setw(6) + << (k == 0 ? 0 : last_l2_error / Velocity_L2_error) << std::endl + << " H1 error: " << std::setw(12) << Velocity_H1_error + << std::setw(0) << " H1_Conv_rate: " << std::setw(6) + << (k == 0 ? 0 : last_H1_error / Velocity_H1_error) << std::endl + << " Hdiv error1: " << std::setw(12) << Velocity_Hdiv_error1 + << std::setw(0) << " Hdiv_Conv_rate1: " << std::setw(6) + << (k == 0 ? 0 : last_Hdiv_error1 / Velocity_Hdiv_error1) + << std::endl + << " L2 pressure: " << std::setw(12) << Pressure_L2_error + << std::setw(0) << " rate: " << std::setw(6) + << (k == 0 ? 0 : last_Pressure_L2_error / Pressure_L2_error) + << std::endl + << std::setw(0) << std::endl; + last_l2_error = Velocity_L2_error; + last_H1_error = Velocity_H1_error; + last_Hdiv_error1 = Velocity_Hdiv_error1; + last_Pressure_L2_error = Pressure_L2_error; + } + + + + // @sect4{StokesProblem::run} + + // The last step in the Stokes class is, as usual, the function that + // generates the initial grid and calls the other functions in the + // respective order. + template + void + StokesProblem::run() + { + GridGenerator::hyper_cube(triangulation); + + triangulation.refine_global(1); + + deallog << " Now running with " << fe.get_name() << std::endl; + + for (unsigned int refinement_cycle = 0; refinement_cycle < 5; + ++refinement_cycle) + { + if (refinement_cycle > 0) + triangulation.refine_global(); + + setup_dofs(); + + int assemble_type = 2; + + assemble_system_mesh_loop(); + + solve(); + + compute_errors(refinement_cycle); + } + } +} // namespace StokesTests + + + +int +main() +{ + try + { + using namespace dealii; + using namespace StokesTests; + const int dim = 2; + + initlog(); + + std::unique_ptr> fe; + const int degree = 2; + + Assert(degree >= 1, ExcMessage("invalid degree!")); + fe = std::make_unique>( + FESystem(FE_DGQ(degree), dim), 1, FE_DGQ(degree - 1), 1); + + deallog << fe->get_name() << ": degree=" << fe->degree + << " tensor_degree=" << fe->tensor_degree() << std::endl; + StokesProblem flow_problem(*fe.get(), degree); + + flow_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; + } + + return 0; +} diff --git a/tests/feinterface/stokes.output b/tests/feinterface/stokes.output new file mode 100644 index 0000000000..6f6be18c16 --- /dev/null +++ b/tests/feinterface/stokes.output @@ -0,0 +1,43 @@ + +DEAL::FESystem<2>[FESystem<2>[FE_DGQ<2>(2)^2]-FE_DGQ<2>(1)]: degree=2 tensor_degree=2 +DEAL:: Now running with FESystem<2>[FESystem<2>[FE_DGQ<2>(2)^2]-FE_DGQ<2>(1)] +DEAL:: Number of active cells: 4 +DEAL:: Number of degrees of freedom: 88 (72+16) +DEAL:: At 1th mesh +DEAL:: L2 error: 0.194224 L2_Conv_rate: 0.00000 +DEAL:: H1 error: 3.65712 H1_Conv_rate: 0.00000 +DEAL:: Hdiv error1: 2.93726 Hdiv_Conv_rate1: 0.00000 +DEAL:: L2 pressure: 0.306514 rate: 0.00000 +DEAL:: +DEAL:: Number of active cells: 16 +DEAL:: Number of degrees of freedom: 352 (288+64) +DEAL:: At 2th mesh +DEAL:: L2 error: 0.0397041 L2_Conv_rate: 4.89179 +DEAL:: H1 error: 1.27499 H1_Conv_rate: 2.86835 +DEAL:: Hdiv error1: 0.624393 Hdiv_Conv_rate1: 4.70418 +DEAL:: L2 pressure: 0.116859 rate: 2.62293 +DEAL:: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 1408 (1152+256) +DEAL:: At 3th mesh +DEAL:: L2 error: 0.00512759 L2_Conv_rate: 7.74323 +DEAL:: H1 error: 0.322333 H1_Conv_rate: 3.95551 +DEAL:: Hdiv error1: 0.159620 Hdiv_Conv_rate1: 3.91175 +DEAL:: L2 pressure: 0.0134970 rate: 8.65819 +DEAL:: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 5632 (4608+1024) +DEAL:: At 4th mesh +DEAL:: L2 error: 0.000650495 L2_Conv_rate: 7.88259 +DEAL:: H1 error: 0.0806845 H1_Conv_rate: 3.99499 +DEAL:: Hdiv error1: 0.0400590 Hdiv_Conv_rate1: 3.98461 +DEAL:: L2 pressure: 0.00186869 rate: 7.22268 +DEAL:: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 22528 (18432+4096) +DEAL:: At 5th mesh +DEAL:: L2 error: 8.19268e-05 L2_Conv_rate: 7.93996 +DEAL:: H1 error: 0.0201660 H1_Conv_rate: 4.00102 +DEAL:: Hdiv error1: 0.0100240 Hdiv_Conv_rate1: 3.99631 +DEAL:: L2 pressure: 0.000321058 rate: 5.82042 +DEAL::