From 99968e929a2e3e4bdd5e30ada1d50f9f8cf44bbe Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 8 Jul 2011 01:36:31 +0000 Subject: [PATCH] Replace a quadratic algorithm by a linear one. This greatly accelerates a good number of the tests we have. For example, deal.II/grid_in and deal.II/grid_in_msh_version_2 speed up from several minutes to just a few seconds. git-svn-id: https://svn.dealii.org/trunk@23935 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 ++ .../deal.II/grid/grid_reordering_internal.h | 18 +-- deal.II/source/grid/grid_reordering.cc | 118 +++++++++++------- 3 files changed, 82 insertions(+), 61 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1dbab530f8..049f03b4f4 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -216,6 +216,13 @@ should be fixed now.

Specific improvements

    +
  1. Fixed: The 2d grid reordering algorithm that is used by all grid readers had +a component that was quadratic in its complexity, sometimes leading to cases +where reading in a mesh in debug mode could take minutes for just a few tens +of thousands of cells. This has now been fixed. +
    +(Wolfgang Bangerth 2011/07/07) +
  2. New: The function DoFTools::count_dofs_per_component now also works for objects of type hp::DoFHandler.
    diff --git a/deal.II/include/deal.II/grid/grid_reordering_internal.h b/deal.II/include/deal.II/grid/grid_reordering_internal.h index c7e5261d19..4f926379a1 100644 --- a/deal.II/include/deal.II/grid/grid_reordering_internal.h +++ b/deal.II/include/deal.II/grid/grid_reordering_internal.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -142,22 +142,6 @@ namespace internal * data of this object. */ CellData<2> original_cell_data; - - /** - * Makes an MQuad from the - * given CellData and MSide - * list. Is derived from - * binary_function to be - * usable with STL - * containers. - * - * Also assumes that the - * edges listed present in - * the CellData are already - * present in the elist - * vector. - */ - struct MakeQuad; }; /** diff --git a/deal.II/source/grid/grid_reordering.cc b/deal.II/source/grid/grid_reordering.cc index 779213bb2b..dbaa9f7190 100644 --- a/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/source/grid/grid_reordering.cc @@ -266,48 +266,74 @@ namespace internal } - struct MQuad::MakeQuad : public std::binary_function, - std::vector, - MQuad> + namespace { - MQuad operator()(const CellData<2> &q, - const std::vector &elist) const - { - // compute the indices of the - // four sides that bound this - // quad - unsigned int edges[4] = { numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int }; - - unsigned int index = 0; - for (std::vector::const_iterator - edge = elist.begin(); edge != elist.end(); ++edge, ++index) - for (unsigned int i=0; i<4; ++i) - // for each of the four - // sides, find the place in - // the list with the - // corresponding line and - // store it. shortcut the - // comparison if it has - // already been found - if ((edges[i] == numbers::invalid_unsigned_int) - && - (MSide::SideSortLess()(*edge, quadside(q,i)) == false)) - edges[i] = index; - - for (unsigned int i=0; i<4; ++i) - Assert (edges[i] != numbers::invalid_unsigned_int, - ExcInternalError()); + /** + * A replacement for std::lower_bound + * if we know that the input sequence + * is in fact sorted. In that case, we + * can do a binary search. + */ + template + Iterator + lower_bound_in_sorted_array (Iterator first, Iterator last, + const T& val, Comp comp) + { + unsigned int length = std::distance(first, last); + unsigned int half; + Iterator middle; - return MQuad(q.vertices[0],q.vertices[1], q.vertices[2], q.vertices[3], - edges[0], edges[1], edges[2], edges[3], - q); + while (length > 0) + { + half = length >> 1; + middle = first; + std::advance(middle, half); + if (comp(*middle, val)) + { + first = middle; + ++first; + length = length - half - 1; + } + else + length = half; } + return first; + } + - }; - + /** + * Create an MQuad object from the + * indices of the four vertices by + * looking up the indices of the four + * sides. + */ + MQuad build_quad_from_vertices(const CellData<2> &q, + const std::vector &elist) + { + // compute the indices of the four + // sides that bound this quad. note + // that the incoming list elist is + // sorted with regard to the + // MSide::SideSortLess criterion + unsigned int edges[4] = { numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int }; + + for (unsigned int i=0; i<4; ++i) + edges[i] = (lower_bound_in_sorted_array (elist.begin(), + elist.end(), + quadside(q,i), + MSide::SideSortLess()) + - + elist.begin()); + + return MQuad(q.vertices[0],q.vertices[1], q.vertices[2], q.vertices[3], + edges[0], edges[1], edges[2], edges[3], + q); + } + } + void @@ -324,7 +350,6 @@ namespace internal { //Reserve some space sides.reserve(4*inquads.size()); - mquads.reserve(inquads.size()); //Insert all the sides into the side vector for (int i = 0;i<4;++i) @@ -348,11 +373,16 @@ namespace internal // Swap trick to shrink the // side vector std::vector(sides).swap(sides); - - //Assigns the correct sides to - //each quads - std::transform(inquads.begin(),inquads.end(), std::back_inserter(mquads), - std::bind2nd(MQuad::MakeQuad(),sides) ); + + // Now assign the correct sides to + // each quads + mquads.reserve(inquads.size()); + std::transform(inquads.begin(), + inquads.end(), + std::back_inserter(mquads), + std_cxx1x::bind(build_quad_from_vertices, + std_cxx1x::_1, + std_cxx1x::cref(sides)) ); // Assign the quads to their sides also. int qctr = 0; -- 2.39.5