From: Guido Kanschat Date: Fri, 23 Aug 2013 20:20:53 +0000 (+0000) Subject: clean code in Assembler::MGMatrixSimple X-Git-Tag: v8.1.0~980 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47eb23ba26c99fb29d390fa18e2c2402a63707f9;p=dealii.git clean code in Assembler::MGMatrixSimple add mg::SparseMatrixCollection git-svn-id: https://svn.dealii.org/trunk@30457 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h index e7a9c73544..e43257df15 100644 --- a/deal.II/include/deal.II/meshworker/simple.h +++ b/deal.II/include/deal.II/meshworker/simple.h @@ -256,29 +256,24 @@ namespace MeshWorker { public: /** - * Constructor, initializing - * the #threshold, which - * limits how small numbers - * may be to be entered into - * the matrix. + * Constructor, initializing the #threshold, which limits how + * small numbers may be to be entered into the matrix. */ MGMatrixSimple(double threshold = 1.e-12); /** - * Store the result matrix - * for later assembling. + * Store the result matrix for later assembling. */ void initialize(MGLevelObject &m); /** - * Initialize the multilevel - * constraints. + * Initialize the multilevel constraints. */ void initialize(const MGConstrainedDoFs &mg_constrained_dofs); /** - * @deprecated This function is of no effect. Only the block info - * structure in DoFInfo is being used. + * @deprecated This function is of no effect. Only the block + * info structure in DoFInfo is being used. * * Store information on the local block structure. If the * assembler is inititialized with this function, @@ -920,8 +915,7 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - if (//mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + if (!mg_constrained_dofs->at_refinement_edge(level, i2[k])) G.add(i1[j], i2[k], M(k,j)); } } @@ -950,8 +944,7 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - if (//mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + if (!mg_constrained_dofs->at_refinement_edge(level, i2[k])) G.add(i1[j], i2[k], M(j,k)); } } @@ -967,49 +960,41 @@ namespace MeshWorker { AssertDimension(M.m(), i1.size()); AssertDimension(M.n(), i2.size()); - - if (mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); - } - else - { - for (unsigned int j=0; j= threshold) - // Enter values into matrix only if j corresponds to a - // degree of freedom on the refinemenent edge, k does - // not, and both are not on the boundary. This is part - // the difference between the complete matrix with no - // boundary condition at the refinement edge and and - // the matrix assembled above by assemble(). - - // Thus the logic is: enter the row if it is - // constrained by hanging node constraints (actually, - // the whole refinement edge), but not if it is - // constrained by a boundary constraint. - if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - { - if (mg_constrained_dofs->set_boundary_values()) - { - if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && - i1[j] == i2[k])) - G.add(i1[j], i2[k], M(j,k)); - } - else - G.add(i1[j], i2[k], M(j,k)); - } - } + Assert(mg_constrained_dofs != 0, ExcInternalError()); + + for (unsigned int j=0; j= threshold) + // Enter values into matrix only if j corresponds to a + // degree of freedom on the refinemenent edge, k does + // not, and both are not on the boundary. This is part + // the difference between the complete matrix with no + // boundary condition at the refinement edge and and + // the matrix assembled above by assemble(). + + // Thus the logic is: enter the row if it is + // constrained by hanging node constraints (actually, + // the whole refinement edge), but not if it is + // constrained by a boundary constraint. + if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + else + G.add(i1[j], i2[k], M(j,k)); + } } - + + template inline void MGMatrixSimple::assemble_out( @@ -1021,38 +1006,29 @@ namespace MeshWorker { AssertDimension(M.n(), i1.size()); AssertDimension(M.m(), i2.size()); - - if (mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(k,j)); - } - else - { - for (unsigned int j=0; j= threshold) - if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - { - if (mg_constrained_dofs->set_boundary_values()) - { - if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && - i1[j] == i2[k])) - G.add(i1[j], i2[k], M(k,j)); - } - else - G.add(i1[j], i2[k], M(k,j)); - } - } + Assert(mg_constrained_dofs != 0, ExcInternalError()); + + for (unsigned int j=0; j= threshold) + if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(k,j)); + } + else + G.add(i1[j], i2[k], M(k,j)); + } } - + template template diff --git a/deal.II/include/deal.II/multigrid/sparse_matrix_collection.h b/deal.II/include/deal.II/multigrid/sparse_matrix_collection.h new file mode 100644 index 0000000000..a7dd03816e --- /dev/null +++ b/deal.II/include/deal.II/multigrid/sparse_matrix_collection.h @@ -0,0 +1,124 @@ +// --------------------------------------------------------------------- +// $Id: mg_matrix.h 30036 2013-07-18 16:55:32Z maier $ +// +// Copyright (C) 2003 - 2013 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. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__mg_sparse_matrix_collection_h +#define __deal2__mg_sparse_matrix_collection_h + +#include +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace mg +{ +/** + * Handler and storage for all five SparseMatrix object involved in + * using multigrid with local refinement. + * + * @author Baerbel Janssen, Guido Kanschat + * @date 2013 + */ + template + class SparseMatrixCollection : public Subscriptor + { + public: + void resize(const unsigned int minlevel, const unsigned int maxlevel); + + template + void reinit(const DH& dof_handler); + + void set_zero(); + + MGLevelObject sparsity; + MGLevelObject sparsity_edge; + + MGLevelObject > matrix; + MGLevelObject > matrix_down; + MGLevelObject > matrix_up; + MGLevelObject > matrix_in; + MGLevelObject > matrix_out; + }; + + + template + void + SparseMatrixCollection::resize(const unsigned int minlevel, const unsigned int maxlevel) + { + matrix.resize(minlevel, maxlevel); + matrix.clear(); + matrix_up.resize(minlevel+1, maxlevel); + matrix_up.clear(); + matrix_down.resize(minlevel+1, maxlevel); + matrix_down.clear(); + matrix_in.resize(minlevel, maxlevel); + matrix_in.clear(); + matrix_out.resize(minlevel, maxlevel); + matrix_out.clear(); + sparsity.resize(minlevel, maxlevel); + sparsity_edge.resize(minlevel, maxlevel); + } + + + template + template + void + SparseMatrixCollection::reinit(const DH& dof_handler) + { + AssertIndexRange(sparsity.max_level(), dof_handler.get_tria().n_levels()); + + for (unsigned int level=sparsity.min_level(); + level<=sparsity.max_level();++level) + { + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level); + sparsity[level].copy_from(c_sparsity); + matrix[level].reinit(sparsity[level]); + matrix_in[level].reinit(sparsity[level]); + matrix_out[level].reinit(sparsity[level]); + if (level>0) + { + CompressedSparsityPattern ci_sparsity; + ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level); + sparsity_edge[level].copy_from(ci_sparsity); + matrix_up[level].reinit(sparsity_edge[level]); + matrix_down[level].reinit(sparsity_edge[level]); + } + } + } + + template + void + SparseMatrixCollection::set_zero() + { + matrix = 0.; + matrix_in = 0.; + matrix_out = 0.; + matrix_up = 0.; + matrix_down = 0.; + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif