From 300afefa416daafed8fd1e91313b9a779047cf09 Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 7 May 2014 17:17:14 +0000 Subject: [PATCH] Added QSorted quadrature rule. git-svn-id: https://svn.dealii.org/trunk@32890 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 8 + deal.II/doc/users/config.sample | 83 +++--- deal.II/include/deal.II/base/quadrature_lib.h | 23 ++ deal.II/source/base/quadrature_lib.cc | 29 +++ tests/base/quadrature_sorted_test.cc | 246 ++++++++++++++++++ tests/base/quadrature_sorted_test.output | 148 +++++++++++ 6 files changed, 486 insertions(+), 51 deletions(-) create mode 100644 tests/base/quadrature_sorted_test.cc create mode 100644 tests/base/quadrature_sorted_test.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index a6597cf1d4..e52d5cf9ad 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,14 @@ inconvenience this causes.

Specific improvements

    +
  1. New: There is now a QSorted quadrature which takes an + arbitrary quadrature at construction time and reorders the quadrature + points according to the weigths, from smaller to bigger. This should + improve stability of higher order polynomial integration. +
    + (Luca Heltai, 2014/05/07) +
    2. +
  2. New: The class VectorizedArray now provides methods VectorizedArray::load(ptr) to read from arbitrary pointer addresses and VectorizedArray::store(ptr) to write to arbitrary pointer addresses, diff --git a/deal.II/doc/users/config.sample b/deal.II/doc/users/config.sample index a64c3dcbb4..b37132c370 100644 --- a/deal.II/doc/users/config.sample +++ b/deal.II/doc/users/config.sample @@ -1,37 +1,9 @@ -## --------------------------------------------------------------------- -## $Id$ -## -## Copyright (C) 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. -## -## --------------------------------------------------------------------- - -# -# Example configuration file -# -# This file can be used as a script to preseed the initial cache of a build -# directory. To do so, invoke cmake with -# -# $ cmake -C configuration.cmake [...] -# -# -# Please note that this file is actually intrepeted as a full featured -# CMake script file, so in order to set variables in the cache a statement -# of the form -# -# SET( CACHE "") -# -# has to be used. See doc/readme.html and doc/users/cmake.html for further -# details on how to use the cmake build system of deal.II. -# +## ## +# Example configuration file # +# # +# See doc/readme.html and doc/users/cmake.html for further # +# details on how to use the cmake build system of deal.II. # +## ## ########################################################################### @@ -119,8 +91,8 @@ # "Fortran Compiler." # ) # -# SET(DEAL_II_CXX_FLAGS "" CACHE STRING -# "The user supplied cache variable will be appended _at the end_ of the auto generated DEAL_II_CXX_FLAGS variable" +# SET(CMAKE_CXX_FLAGS "" CACHE STRING +# "The user supplied cache variable will be appended _at the end_ of the auto generated CMAKE_CXX_FLAGS variable" # ) # # SET(DEAL_II_CXX_FLAGS_DEBUG "" CACHE STRING @@ -143,18 +115,6 @@ # "The user supplied cache variable will be appended _at the end_ of the auto generated DEAL_II_LINKER_FLAGS_RELEASE variable" # ) # -# SET(DEAL_II_DEFINITIONS "" CACHE STRING -# "Additional, user supplied compile definitions" -# ) -# -# SET(DEAL_II_DEFINITIONS_DEBUG "" CACHE STRING -# "Additional, user supplied compile definitions" -# ) -# -# SET(DEAL_II_DEFINITIONS_RELEASE "" CACHE STRING -# "Additional, user supplied compile definitions" -# ) -# # SET(BUILD_SHARED_LIBS "ON" CACHE BOOL # "Build a shared library" # ) @@ -175,6 +135,7 @@ # "If set to ON, then use 64-bit data types to represent global degree of freedom indices. The default is to OFF. You only want to set this to ON if you will solve problems with more than 2^31 (approximately 2 billion) unknowns. If set to ON, you also need to ensure that both Trilinos and/or PETSc support 64-bit indices." # ) # +# ########################################################################### @@ -415,6 +376,25 @@ # +# +# muPaser: +# +# SET(DEAL_II_WITH_MUPARSER ON CACHE BOOL +# "Build deal.II with support for muparser" +# +# Automatic detection: +# +# Specify a hint with CMAKE_PREFIX_PATH or by setting +# SET(MUPARSER_DIR "/.../..." CACHE PATH "") +# +# Manual setup: +# +# SET(MUPARSER_FOUND TRUE CACHE BOOL "") +# SET(MUPARSER_LIBRARIES "library;and;semicolon;separated;list;of;link;interface" CACHE STRING "") +# SET(MUPARSER_INCLUDE_DIRS "semicolon;separated;list;of;include;dirs" CACHE STRING "") +# + + # # Netcdf: # @@ -632,11 +612,12 @@ # C++11 support is autodetected. You can explicitly disable C+11 support by # specifying # -# SET(DEAL_II_WITH_CXX11 FALSE CACHE BOOL "") +# SET(DEAL_II_HAVE_CXX11_FLAG FALSE CACHE BOOL "") # # A custom C++11 flag can be set by setting # -# SET(DEAL_II_CXX11_FLAG "-std=c++11" CACHE STRING "") +# SET(DEAL_II_HAVE_CXX11_FLAG TRUE CACHE BOOL "") +# SET(DEAL_II_CXX11_FLAG "-std=c++0x" CACHE STRING "") # @@ -674,7 +655,7 @@ # "Library suffix for the debug library" # ) # -# SET(DEAL_II_RELEASE_SUFFIX "" CACHE STRING +# SET_IF_EMPTY(DEAL_II_RELEASE_SUFFIX "" CACHE STRING # "Library suffix for the release library" # ) # diff --git a/deal.II/include/deal.II/base/quadrature_lib.h b/deal.II/include/deal.II/base/quadrature_lib.h index a448a80c9e..5bfe461aa0 100644 --- a/deal.II/include/deal.II/base/quadrature_lib.h +++ b/deal.II/include/deal.II/base/quadrature_lib.h @@ -435,6 +435,29 @@ private: +/** + * Sorted Quadrature. Given an arbitrary quadrature formula, this + * class generates a quadrature formula where the quadrature points + * are ordered according the weights, from those with smaller + * corresponding weight, to those with higher corresponding weights. + * This might be necessary, for example, when integrating high order + * polynomials, since in these cases you might sum very big numbers + * with very small numbers, and summation is not stable if the numbers + * to sum are not close to each other. + */ +template +class QSorted : public Quadrature +{ +public: + /** The constructor takes an arbitrary quadrature formula. */ + QSorted (const Quadrature); + + /** A rule to reorder pairs of points and weights.*/ + bool operator()(const std::pair > &a, + const std::pair > &b); +}; + + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ diff --git a/deal.II/source/base/quadrature_lib.cc b/deal.II/source/base/quadrature_lib.cc index 7752b9d81e..8bebce9102 100644 --- a/deal.II/source/base/quadrature_lib.cc +++ b/deal.II/source/base/quadrature_lib.cc @@ -936,6 +936,31 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, } +template +QSorted::QSorted(Quadrature quad) : + Quadrature(quad.size()) +{ + std::vector< std::pair > > wp; + for (unsigned int i=0; i >(quad.weight(i), + quad.point(i))); + sort(wp.begin(), wp.end(), *this); + for(unsigned int i=0; iweights[i] = wp[i].first; + this->quadrature_points[i] = wp[i].second; + } +} + + +template +bool QSorted::operator()(const std::pair > &a, + const std::pair > &b) +{ + return (a.first < b.first); +} + + // construct the quadrature formulae in higher dimensions by // tensor product of lower dimensions @@ -1009,4 +1034,8 @@ template class QSimpson<3>; template class QMilne<3>; template class QWeddle<3>; +template class QSorted<1>; +template class QSorted<2>; +template class QSorted<3>; + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/quadrature_sorted_test.cc b/tests/base/quadrature_sorted_test.cc new file mode 100644 index 0000000000..6393b765e1 --- /dev/null +++ b/tests/base/quadrature_sorted_test.cc @@ -0,0 +1,246 @@ +// --------------------------------------------------------------------- +// $Id: quadrature_test.cc 31349 2013-10-20 19:07:06Z maier $ +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + + + +// sort quadratures according to weights, and verify that they +// integrate accurately polynomial of increasing order. +// Verify that weights are actually sorted + + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include + +template +void +fill_vector (std::vector *> &quadratures) +{ + quadratures.push_back (new QSorted(QMidpoint())); + quadratures.push_back (new QSorted(QTrapez())); + quadratures.push_back (new QSorted(QSimpson())); + quadratures.push_back (new QSorted(QMilne())); + quadratures.push_back (new QSorted(QWeddle())); + for (unsigned int i=0; i<9; ++i) + { + quadratures.push_back (new QSorted(QGauss(i))); + } + QMilne<1> q1d; + quadratures.push_back (new QSorted(Quadrature(q1d))); + for (unsigned int i=2; i<8; ++i) + { + quadratures.push_back (new QSorted(QGaussLobatto(i))); + } +} + +template +void +check_cells (std::vector*> &quadratures) +{ + Quadrature quadrature; + for (unsigned int n=0; n > &points=quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + deallog << "Quadrature no." << n; + + unsigned int i=0; + double quadrature_int=0; + double exact_int=0; + double err = 0; + + do + { + ++i; + + quadrature_int=0; + // Check the polynomial x^i*y^i + + for (unsigned int x=0; x(points[x](2)), i*1.0); + case 2: + f *= std::pow(static_cast(points[x](1)), i*1.0); + case 1: + f *= std::pow(static_cast(points[x](0)), i*1.0); + } + quadrature_int+=f*weights[x]; + } + + // the exact integral is 1/(i+1) + exact_int=1./std::pow(static_cast(i+1),dim); + err = std::fabs(quadrature_int-exact_int); + } + while (err<1e-14); + // Uncomment here for testing +// deallog << " (Int " << quadrature_int << ',' << exact_int << ")"; + deallog << " is exact for polynomials of degree " << i-1 << std::endl; + } +} + + +template +void +check_faces (const std::vector*>& quadratures, const bool sub) +{ + if (sub) + deallog.push("subfaces"); + else + deallog.push("faces"); + + for (unsigned int n=0; n quadrature (sub == false? + QProjector::project_to_all_faces(*quadratures[n]) : + QProjector::project_to_all_subfaces(*quadratures[n])); + const std::vector > &points=quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + deallog << "Quadrature no." << n; + + unsigned int i=0; + long double quadrature_int=0; + double exact_int=0; + double err = 0; + + do + { + ++i; + + quadrature_int=0; + // Check the polynomial + // x^i*y^i*z^i + + for (unsigned int x=0; x