From: ESeNonFossiIo Date: Mon, 7 Sep 2015 14:37:55 +0000 (+0200) Subject: constraint linear operators X-Git-Tag: v8.4.0-rc2~126^2~12 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=abd3b2799dd4773df1080df74e562f47e46689f3;p=dealii.git constraint linear operators --- diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index eccf331ab1..8533d91186 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -1116,6 +1116,98 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix) return return_op; } +template , + typename Domain = Range, + typename Matrix> +LinearOperator +constraints_linear_operator(const ConstraintMatrix &, const Matrix &); + +template , + typename Domain = Range, + typename Matrix> +LinearOperator +constrained_linear_operator(const ConstraintMatrix &, const Matrix &); + +/** + * @relates LinearOperator + * + * Description... TODO! + * + * @ingroup LAOperators + */ +template < typename Range, + typename Domain, + typename Matrix > +LinearOperator +constraints_linear_operator(const ConstraintMatrix &cm, const Matrix &m) +{ + LinearOperator return_op = linear_operator(m); + + return_op.vmult = [&cm](Range &v, const Domain &u) + { + for (auto i : v.locally_owned_elements()) + { + if (cm.is_constrained(i)) + { + v(i) = 0; + const std::vector< std::pair < types::global_dof_index, double > > + *entries = cm.get_constraint_entries (i); + for (types::global_dof_index j=0; j < entries->size(); ++j) + { + types::global_dof_index pos = (*entries)[j].first; + v(i) += u(pos) * (*entries)[j].second; + } + } + else + v(i) = u(i); + } + }; + + return_op.Tvmult = [&cm](Range &v, const Domain &u) + { + for (auto i : u.locally_owned_elements()) + { + if (cm.is_constrained(i)) + { + v(i)=0; + const std::vector< std::pair < types::global_dof_index, double > > + *entries = cm.get_constraint_entries (i); + for (types::global_dof_index j=0; j < entries->size(); ++j) + { + types::global_dof_index pos = (*entries)[j].first; + v(pos) += u(i) * (*entries)[j].second; + } + } + else + v(i)=u(i); + + } + }; + + return return_op; +} + +/** + * @relates LinearOperator + * + * Description... TODO! + * + * @ingroup LAOperators + */ +template < typename Range, + typename Domain, + typename Matrix > +LinearOperator +constrained_linear_operator(const ConstraintMatrix &cm, const Matrix &m) +{ + auto C = constraints_linear_operator(cm, m); + auto Ct = transpose_operator(C); + auto A = linear_operator(m); + + return Ct * A * C; +} + + //@} DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/packaged_operation.h b/include/deal.II/lac/packaged_operation.h index 76d3144ffe..933eb1dfe8 100644 --- a/include/deal.II/lac/packaged_operation.h +++ b/include/deal.II/lac/packaged_operation.h @@ -827,6 +827,48 @@ operator*(const PackagedOperation &comp, return return_comp; } +template , + typename Domain = Range, + typename Matrix> +PackagedOperation +constrained_rhs(const ConstraintMatrix &, const Matrix &, const Range &); + +/** + * @relates LinearOperator + * + * Description... TODO! + * + * @ingroup LAOperators + */ +// template class PackagedOperation; +template +PackagedOperation +constrained_rhs(const ConstraintMatrix &cm, const Matrix &m, const Range &rhs) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = [&m](Range &v, bool fast) + { + internal::LinearOperator::ReinitHelper::reinit_range_vector(m, v, fast); + }; + + return_comp.apply = [&cm, &m, &rhs](Range &v) + { + const auto M = linear_operator(m); + const auto C = constraints_linear_operator(cm, m); + const auto Ct = transpose_operator(C); + + for (auto i : v.locally_owned_elements()) + v[i] = -1 * cm.get_inhomogeneity(i); + + v = Ct * ( rhs + M * v); + }; + + return return_comp; +} + //@} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/linear_operator_09.cc b/tests/lac/linear_operator_09.cc new file mode 100644 index 0000000000..6369b81617 --- /dev/null +++ b/tests/lac/linear_operator_09.cc @@ -0,0 +1,482 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 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. +// +// --------------------------------------------------------------------- + +// This test is based on step-6 +// The aim is to test constrained linear operators +// and compare them to the old implementation. + +#include "../tests.h" + +#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 + +using namespace dealii; + +template +class Step6 +{ +public: + Step6 (); + ~Step6 (); + + void run (); + +private: + void setup_system (); + + void assemble_system (); + void assemble_system_lo (); + + void solve (); + void solve_lo (); + + void refine_grid (); + + Triangulation triangulation; + Triangulation triangulation_lo; + + DoFHandler dof_handler; + DoFHandler dof_handler_lo; + + FE_Q fe; + FE_Q fe_lo; + + ConstraintMatrix constraints; + ConstraintMatrix constraints_lo; + + SparsityPattern sparsity_pattern; + SparsityPattern sparsity_pattern_lo; + + SparseMatrix system_matrix; + SparseMatrix system_matrix_lo; + + Vector solution; + Vector solution_lo; + + Vector system_rhs; + Vector system_rhs_lo; +}; + + + + +template +class Coefficient : public Function +{ +public: + Coefficient () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i +Step6::Step6 () + : + dof_handler (triangulation), + dof_handler_lo (triangulation_lo), + fe (2), + fe_lo (2) +{} + +template +Step6::~Step6 () +{ + dof_handler.clear (); +} + + + +template +void Step6::setup_system () +{ + dof_handler.distribute_dofs (fe); + + solution.reinit (dof_handler.n_dofs()); + solution_lo.reinit (dof_handler.n_dofs()); + + system_rhs.reinit (dof_handler.n_dofs()); + system_rhs_lo.reinit (dof_handler.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + + + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ConstantFunction(1.), + constraints); + + + constraints.close (); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /*keep_constrained_dofs = */ false); + + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit (sparsity_pattern); + + constraints_lo.clear (); + constraints_lo.close (); + + DynamicSparsityPattern dsp_lo(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp_lo, + constraints_lo, + /*keep_constrained_dofs = */ false); + + sparsity_pattern_lo.copy_from(dsp_lo); + + system_matrix_lo.reinit (sparsity_pattern_lo); +} + + + +template +void Step6::assemble_system () +{ + const QGauss quadrature_formula(3); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } +} + +template +void Step6::assemble_system_lo () +{ + const QGauss quadrature_formula(3); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); + constraints_lo.distribute_local_to_global (cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix_lo, + system_rhs_lo); + } + +} + + +template +void Step6::solve () +{ + // constraints.condense(system_matrix, system_rhs); + unsigned int n_iterations = 0; + auto M = linear_operator(system_matrix); + + SolverControl solver_control (10000, 1e-12); + SolverCG<> solver (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + auto M_inv = inverse_operator( M, + solver, + preconditioner); + + + + M_inv.vmult(solution, system_rhs); + n_iterations = solver_control.last_step(); + constraints.distribute (solution); +} + + +template +void Step6::solve_lo () +{ + unsigned int n_iterations = 0; + + auto new_system_rhs_lo = constrained_rhs< >( + constraints, system_matrix_lo, system_rhs_lo);; + auto M = constrained_linear_operator< >( constraints, system_matrix_lo ); + + SolverControl solver_control (10000, 1e-12); + SolverCG<> solver (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + auto M_inv = inverse_operator( M, + solver, + preconditioner); + + M_inv.vmult(solution_lo, new_system_rhs_lo); + n_iterations = solver_control.last_step(); + + constraints.distribute(solution_lo); +} + + + +template +void Step6::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + triangulation.execute_coarsening_and_refinement (); +} + +template +void Step6::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_ball (triangulation); + + static const SphericalManifold boundary; + triangulation.set_all_manifold_ids_on_boundary(0); + triangulation.set_manifold (0, boundary); + + triangulation.refine_global (1); + } + else + refine_grid (); + + setup_system (); + + assemble_system (); + solve (); + + assemble_system_lo (); + solve_lo (); + + // auto rhs_tmp = system_rhs; + // new_system_rhs_lo.apply(rhs_tmp); + // compare_solution(system_rhs, rhs_tmp, dof_handler.n_dofs(), 1e-6, false, false, "RHS"); + Vector diff = solution-solution_lo; + if (diff.l2_norm() < 1e-10) + deallog << "OK" << std::endl; + } +} + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + { + try + { + Step6<2> laplace_problem_2d; + laplace_problem_2d.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; + + deallog << "OK" << std::endl; + } +} diff --git a/tests/lac/linear_operator_09.with_cxx11=on.output b/tests/lac/linear_operator_09.with_cxx11=on.output new file mode 100644 index 0000000000..80feeb7ee1 --- /dev/null +++ b/tests/lac/linear_operator_09.with_cxx11=on.output @@ -0,0 +1,16 @@ + +DEAL:cg::Starting value 13.25516422 +DEAL:cg::Convergence step 24 value 1.400407060e-13 +DEAL:cg::Starting value 13.25516422 +DEAL:cg::Convergence step 24 value 1.399854635e-13 +DEAL::OK +DEAL:cg::Starting value 17.40605136 +DEAL:cg::Convergence step 33 value 3.786179069e-13 +DEAL:cg::Starting value 17.40605136 +DEAL:cg::Convergence step 33 value 3.789771041e-13 +DEAL::OK +DEAL:cg::Starting value 17.40587637 +DEAL:cg::Convergence step 56 value 4.293614863e-13 +DEAL:cg::Starting value 17.40587637 +DEAL:cg::Convergence step 56 value 4.295516840e-13 +DEAL::OK