From: Timo Heister Date: Fri, 5 Jul 2024 10:13:20 +0000 (-0400) Subject: add tests/feinterface/stokes_simplex X-Git-Tag: v9.6.0-rc1~120^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=57e6e6fc5db72dcfc12f39b135dd69dcce704853;p=dealii.git add tests/feinterface/stokes_simplex This is adding a meshworker simplex test for DG Stokes, extracted from #13868 --- diff --git a/tests/feinterface/stokes_simplex.cc b/tests/feinterface/stokes_simplex.cc new file mode 100644 index 0000000000..ae3189cff4 --- /dev/null +++ b/tests/feinterface/stokes_simplex.cc @@ -0,0 +1,923 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2019 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// test FEInterfaceValues for a DG Stokes problem on a +// simplex mesh. + +#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 +#include +#include + +#include +#include + +#include "../tests.h" + +#define USE_SIMPLEX + +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)); + } + } + + 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; + +#ifdef USE_SIMPLEX + const QGaussSimplex quadrature; + const QGaussSimplex face_quadrature; + const QGaussSimplex quadrature_2; + const QGaussSimplex face_quadrature_2; + const MappingFE mapping; +#else + const QGauss quadrature; + const QGauss face_quadrature; + const QGauss quadrature_2; + const QGauss face_quadrature_2; + const MappingQ1 mapping; +#endif + + Triangulation triangulation; + FiniteElement &fe; + DoFHandler dof_handler; + + AffineConstraints<> constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_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) + , quadrature(pressure_degree + 1) + , face_quadrature(pressure_degree + 1) + , quadrature_2(pressure_degree + 2) + , face_quadrature_2(pressure_degree + 2) +#ifdef USE_SIMPLEX + , mapping(FE_SimplexP(1)) +#else + , mapping() +#endif + , fe(fe) + , dof_handler(triangulation) + {} + + template + void + StokesProblem::setup_dofs() + { + system_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); + DoFRenumbering::block_wise(dof_handler); + + std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, 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 ©_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)) * + 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 ©_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->measure() / cell->face(face_no)->measure(); + 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 ©_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 ©_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_of_gradients(j, point) * + normals[point]) * + fe_fv[velocities].jump_in_values(i, point) + + // - nu [u] . {\nabla v}n (symmetry) // NIPG: use + + - nu * fe_fv[velocities].jump_in_values(j, point) * + (fe_fv[velocities].average_of_gradients(i, point) * + normals[point]) + + // nu sigma [u].[v] (penalty) + + nu * penalty * + scalar_product(fe_fv[velocities].jump_in_values(j, point), + fe_fv[velocities].jump_in_values(i, point)) + + // {p} ([v].n) + + fe_fv[pressure].average_of_values(j, point) * + scalar_product(fe_fv[velocities].jump_in_values(i, point), + normals[point]) + + // {q} ([u].n) + + fe_fv[pressure].average_of_values(i, point) * + scalar_product(fe_fv[velocities].jump_in_values(j, point), + normals[point])) * + JxW[point]; + } + }; + + auto copier = [&](const CopyData &c) { + copy(c, constraints, system_matrix, system_rhs); + }; + + 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; + + + 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); + } + + + + template + void + StokesProblem::solve() + { + solution = 0; + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + + constraints.distribute(solution); + + const double mean_pressure = + VectorTools::compute_mean_value(dof_handler, quadrature_2, solution, dim); + VectorTools::add_constant(solution, dof_handler, dim, -mean_pressure); + } + + + + 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(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + quadrature_2, + VectorTools::L2_norm, + &velocity_mask); + + const double Velocity_L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + quadrature_2, + VectorTools::H1_norm, + &velocity_mask); + + const double Velocity_H1_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Functions::ZeroFunction(dim + 1), + difference_per_cell, + quadrature_2, + VectorTools::Hdiv_seminorm, + &velocity_mask); + + const double Velocity_Hdiv_error1 = difference_per_cell.l2_norm(); + + static double last_Pressure_L2_error = 0; + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + quadrature_2, + 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; + } + + + + template + void + StokesProblem::run() + { +#ifdef USE_SIMPLEX + Triangulation temp; + GridGenerator::hyper_cube(temp); + GridGenerator::convert_hypercube_to_simplex_mesh(temp, triangulation); +#else + GridGenerator::hyper_cube(triangulation); +#endif + + triangulation.refine_global(1); + + deallog << " Now running with " << fe.get_name() << std::endl; + + for (unsigned int refinement_cycle = 0; refinement_cycle < 4; + ++refinement_cycle) + { + if (refinement_cycle > 0) + triangulation.refine_global(); + + setup_dofs(); + + assemble_system_mesh_loop(); + + solve(); + + compute_errors(refinement_cycle); + } + } +} // namespace StokesTests + + +template +void +test_stokes(FiniteElement &fe, const unsigned int pressure_degree) +{ + deallog << fe.get_name() << ": degree=" << fe.degree + << " tensor_degree=" << fe.tensor_degree() << std::endl; + StokesTests::StokesProblem flow_problem(fe, pressure_degree); + flow_problem.run(); +} + +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!")); + { +#ifdef USE_SIMPLEX + fe = std::make_unique>( + FESystem(FE_SimplexDGP(degree), dim), + 1, + FE_SimplexDGP(degree - 1), + 1); +#else + fe = std::make_unique>(FESystem(FE_DGQ(degree), + dim), + 1, + FE_DGQ(degree - 1), + 1); +#endif + test_stokes(*fe, degree); + } + } + 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_simplex.with_umfpack=on.output b/tests/feinterface/stokes_simplex.with_umfpack=on.output new file mode 100644 index 0000000000..54396e33fb --- /dev/null +++ b/tests/feinterface/stokes_simplex.with_umfpack=on.output @@ -0,0 +1,35 @@ + +DEAL::FESystem<2>[FESystem<2>[FE_SimplexDGP<2>(2)^2]-FE_SimplexDGP<2>(1)]: degree=2 tensor_degree=2 +DEAL:: Now running with FESystem<2>[FESystem<2>[FE_SimplexDGP<2>(2)^2]-FE_SimplexDGP<2>(1)] +DEAL:: Number of active cells: 32 +DEAL:: Number of degrees of freedom: 480 (384+96) +DEAL:: At 1th mesh +DEAL:: L2 error: 0.117546 L2_Conv_rate: 0.00000 +DEAL:: H1 error: 2.20951 H1_Conv_rate: 0.00000 +DEAL:: Hdiv error1: 0.430647 Hdiv_Conv_rate1: 0.00000 +DEAL:: L2 pressure: 4.15183 rate: 0.00000 +DEAL:: +DEAL:: Number of active cells: 128 +DEAL:: Number of degrees of freedom: 1920 (1536+384) +DEAL:: At 2th mesh +DEAL:: L2 error: 0.0138896 L2_Conv_rate: 8.46286 +DEAL:: H1 error: 0.577807 H1_Conv_rate: 3.82396 +DEAL:: Hdiv error1: 0.108836 Hdiv_Conv_rate1: 3.95685 +DEAL:: L2 pressure: 1.41599 rate: 2.93210 +DEAL:: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 7680 (6144+1536) +DEAL:: At 3th mesh +DEAL:: L2 error: 0.00130976 L2_Conv_rate: 10.6047 +DEAL:: H1 error: 0.137071 H1_Conv_rate: 4.21539 +DEAL:: Hdiv error1: 0.0288825 Hdiv_Conv_rate1: 3.76823 +DEAL:: L2 pressure: 0.380977 rate: 3.71674 +DEAL:: +DEAL:: Number of active cells: 2048 +DEAL:: Number of degrees of freedom: 30720 (24576+6144) +DEAL:: At 4th mesh +DEAL:: L2 error: 0.000134320 L2_Conv_rate: 9.75105 +DEAL:: H1 error: 0.0329769 H1_Conv_rate: 4.15658 +DEAL:: Hdiv error1: 0.00739000 Hdiv_Conv_rate1: 3.90832 +DEAL:: L2 pressure: 0.0977908 rate: 3.89584 +DEAL::