From: Daniel Arndt Date: Sun, 5 Mar 2017 11:57:41 +0000 (+0100) Subject: changelog entry and a test for a Stokes operator using MatrixFreeOperatorsBase X-Git-Tag: v8.5.0-rc1~69^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=523201f891e3b432339c1afd63c38e4e7d16ca64;p=dealii.git changelog entry and a test for a Stokes operator using MatrixFreeOperatorsBase --- diff --git a/doc/news/changes/minor/20170305DanielArndt b/doc/news/changes/minor/20170305DanielArndt new file mode 100644 index 0000000000..daf28e9009 --- /dev/null +++ b/doc/news/changes/minor/20170305DanielArndt @@ -0,0 +1,4 @@ +Improved: MatrixFreeOperators::Base also supports block vectors. +Derived operators may act on multiple blocks. +
+(Daniel Arndt, 2017/03/05) diff --git a/tests/matrix_free/matrix_vector_stokes_base.cc b/tests/matrix_free/matrix_vector_stokes_base.cc new file mode 100644 index 0000000000..dc68467922 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_base.cc @@ -0,0 +1,292 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + +// Test the correctness of an implementation for the Stokes operator +// based on MatrixFreeOperators::Base by comparing with a matrix version. + +#include "../tests.h" +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MatrixFreeTest : public MatrixFreeOperators::Base +{ +public: + typedef typename BlockVectorType::value_type Number; + typedef typename MatrixFreeOperators::Base Base; + + void compute_diagonal () + { + AssertThrow(false, ExcNotImplemented()); + } + +protected: + + void apply_add(BlockVectorType &dst,const BlockVectorType &src) const + { + Base::data->cell_loop(&MatrixFreeTest::local_apply_cell, this, dst, src); + } + + void local_apply_cell(const dealii::MatrixFree &data, + BlockVectorType &dst, const BlockVectorType &src, + const std::pair &cell_range) const + { + typedef VectorizedArray vector_t; + FEEvaluation velocity (data, 0); + FEEvaluation pressure (data, 1); + + for (unsigned int cell=cell_range.first; cell sym_grad_u = + velocity.get_symmetric_gradient (q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(sym_grad_u); + pressure.submit_value (div, q); + + // subtract p * I + for (unsigned int d=0; d +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_Q fe_u_scal (degree + 1); + FESystem fe_u (fe_u_scal, dim); + FE_Q fe_p(degree); + FESystem fe(fe_u_scal, dim, fe_p, 1); + + DoFHandler dof(tria); + DoFHandler dof_u(tria); + DoFHandler dof_p(tria); + + dof.distribute_dofs(fe); + dof_u.distribute_dofs(fe_u); + dof_p.distribute_dofs(fe_p); + + std::shared_ptr > mf_data; + + dof.distribute_dofs(fe); + ConstraintMatrix constraints, constraints_u, constraints_p; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + LinearAlgebra::distributed::BlockVector mf_system_rhs; + LinearAlgebra::distributed::BlockVector mf_solution; + + std::vector stokes_sub_blocks(dim + 1, 0); + stokes_sub_blocks[dim] = 1; + DoFRenumbering::component_wise(dof, stokes_sub_blocks); + + DoFTools::make_hanging_node_constraints (dof, constraints); + constraints.close(); + DoFTools::make_hanging_node_constraints (dof_u, constraints_u); + constraints_u.close(); + DoFTools::make_hanging_node_constraints (dof_p, constraints_p); + constraints_p.close(); + + std::vector dofs_per_block(2); + DoFTools::count_dofs_per_block(dof, dofs_per_block, stokes_sub_blocks); + { + BlockDynamicSparsityPattern csp(2, 2); + + for (unsigned int d = 0; d < 2; ++d) + { + for (unsigned int e = 0; e < 2; ++e) + csp.block(d, e).reinit(dofs_per_block[d], dofs_per_block[e]); + } + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern(dof, csp, constraints, false); + sparsity_pattern.copy_from(csp); + } + + system_matrix.reinit(sparsity_pattern); + + // this is from step-22 + { + QGauss quadrature_formula(degree + 2); + + FEValues fe_values( + fe, quadrature_formula, update_values | update_JxW_values | update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> phi_grads_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + local_matrix = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + local_matrix(i, j) += (phi_grads_u[i] * phi_grads_u[j] - div_phi_u[i] * phi_p[j] - + phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + } + } + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(local_matrix, local_dof_indices, system_matrix); + } + } + + solution.reinit(2); + mf_solution.reinit(2); + for (unsigned int d = 0; d < 2; ++d) + { + solution.block(d).reinit(dofs_per_block[d]); + mf_solution.block(d).reinit(dofs_per_block[d]); + } + solution.collect_sizes(); + mf_solution.collect_sizes(); + system_rhs.reinit(solution); + mf_system_rhs.reinit(mf_solution); + + // fill system_rhs with random numbers + for (unsigned int j=0; j >(new MatrixFree()); + // setup matrix-free structure + { + std::vector*> dofs; + dofs.push_back(&dof_u); + dofs.push_back(&dof_p); + std::vector constraints; + constraints.push_back (&constraints_u); + constraints.push_back (&constraints_p); + QGauss<1> quad(degree+2); + // no parallelism + mf_data->reinit (dofs, constraints, quad, + typename MatrixFree::AdditionalData + (MPI_COMM_WORLD, + MatrixFree::AdditionalData::none)); + } + system_matrix.vmult (solution, system_rhs); + + MatrixFreeTest > mf; + mf.initialize (mf_data); + mf.vmult (mf_solution, mf_system_rhs); + + // Verification + for (unsigned int i=0; i(); + test<2,2>(); + deallog.pop(); + deallog.push("3D"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + + return 0; +} diff --git a/tests/matrix_free/matrix_vector_stokes_base.output b/tests/matrix_free/matrix_vector_stokes_base.output new file mode 100644 index 0000000000..cc104fc98f --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_base.output @@ -0,0 +1,9 @@ + +DEAL:2D::Verification fe degree 1: 0 +DEAL:2D:: +DEAL:2D::Verification fe degree 2: 0 +DEAL:2D:: +DEAL:3D::Verification fe degree 1: 0 +DEAL:3D:: +DEAL:3D::Verification fe degree 2: 0 +DEAL:3D::