From: Peter Munch Date: Thu, 2 Jul 2020 09:29:29 +0000 (+0200) Subject: Add additional quadrature rule for Simplex X-Git-Tag: v9.3.0-rc1~1341^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10643%2Fhead;p=dealii.git Add additional quadrature rule for Simplex --- diff --git a/doc/news/changes/minor/20200702Munch b/doc/news/changes/minor/20200702Munch new file mode 100644 index 0000000000..b9d4d55d3b --- /dev/null +++ b/doc/news/changes/minor/20200702Munch @@ -0,0 +1,3 @@ +New: A new quadrature rule for simplex geometric entities has been added. +
+(Peter Munch, 2020/07/02) diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h new file mode 100644 index 0000000000..fbe4f88241 --- /dev/null +++ b/include/deal.II/simplex/quadrature_lib.h @@ -0,0 +1,48 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tet_quadrature_lib_h +#define dealii_tet_quadrature_lib_h + + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * Integration rule for simplex entities. + * + * Following number of quadrature points are currently supported: + * - 2D: 1, 3, 7 + * - 3D: 1, 4, 10 + */ + template + class PGauss : public QSimplex + { + public: + /** + * Constructor taking the number of quadrature points @p n_points. + */ + explicit PGauss(const unsigned int n_points); + }; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 7173d25b85..a6d85d8b68 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -71,6 +71,7 @@ ADD_SUBDIRECTORY(differentiation) ADD_SUBDIRECTORY(physics) ADD_SUBDIRECTORY(optimization/rol) ADD_SUBDIRECTORY(non_matching) +ADD_SUBDIRECTORY(simplex) ADD_SUBDIRECTORY(sundials) FOREACH(build ${DEAL_II_BUILD_TYPES}) diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt new file mode 100644 index 0000000000..cfc0183c55 --- /dev/null +++ b/source/simplex/CMakeLists.txt @@ -0,0 +1,42 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2020 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.md at +## the top level directory of deal.II. +## +## --------------------------------------------------------------------- + +INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_unity_include_src + quadrature_lib.cc + ) + +SET(_separate_src + ) + +# include all files in the unity file +SET(_n_includes_per_unity_file 20) + +SETUP_SOURCE_LIST("${_unity_include_src}" + "${_separate_src}" + ${_n_includes_per_unity_file} + _src + ) + +SET(_inst +) + +FILE(GLOB _header + ${CMAKE_SOURCE_DIR}/include/deal.II/tet/*.h + ) + +DEAL_II_ADD_LIBRARY(obj_simplex OBJECT ${_src} ${_header} ${_inst}) +EXPAND_INSTANTIATIONS(obj_simplex "${_inst}") diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc new file mode 100644 index 0000000000..94bd3ad9c6 --- /dev/null +++ b/source/simplex/quadrature_lib.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +#include + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + template + PGauss::PGauss(const unsigned int n_points) + : QSimplex(Quadrature()) + { + // fill quadrature points and quadrature weights + + if (dim == 2) + { + if (n_points == 1) + { + const double p = 1.0 / 3.0; + this->quadrature_points.emplace_back(p, p); + this->weights.emplace_back(1.0); + } + else if (n_points == 3) + { + const double Q23 = 2.0 / 3.0; + const double Q16 = 1.0 / 6.0; + + this->quadrature_points.emplace_back(Q23, Q16); + this->quadrature_points.emplace_back(Q16, Q23); + this->quadrature_points.emplace_back(Q16, Q16); + this->weights.emplace_back(Q16); + this->weights.emplace_back(Q16); + this->weights.emplace_back(Q16); + } + else if (n_points == 7) + { + const double q12 = 0.5; + + // clang-format off + this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330); + this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230); + this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530); + this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230); + this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050); + this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898); + this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050); + // clang-format on + + this->weights.emplace_back(q12 * 0.225); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.132394152789); + this->weights.emplace_back(q12 * 0.132394152789); + this->weights.emplace_back(q12 * 0.132394152789); + } + } + else if (dim == 3) + { + if (n_points == 1) + { + const double Q14 = 1.0 / 4.0; + const double Q16 = 1.0 / 6.0; + + this->quadrature_points.emplace_back(Q14, Q14, Q14); + this->weights.emplace_back(Q16); + } + else if (n_points == 4) + { + const double Q124 = 1.0 / 6.0 / 4.0; + + const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0; + const double pbeta = (5.0 - sqrt(5.0)) / 20.0; + this->quadrature_points.emplace_back(pbeta, pbeta, pbeta); + this->quadrature_points.emplace_back(palpha, pbeta, pbeta); + this->quadrature_points.emplace_back(pbeta, palpha, pbeta); + this->quadrature_points.emplace_back(pbeta, pbeta, palpha); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + } + else if (n_points == 10) + { + const double Q16 = 1.0 / 6.0; + + // clang-format off + this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852); + this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852); + this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444); + this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852); + this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000); + // clang-format on + + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + } + } + + AssertDimension(this->quadrature_points.size(), this->weights.size()); + Assert(this->quadrature_points.size() > 0, + ExcMessage("No valid quadrature points!")); + } + +} // namespace Simplex + + +template class Simplex::PGauss<2>; +template class Simplex::PGauss<3>; + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/simplex/quadrature_lib_01.cc b/tests/simplex/quadrature_lib_01.cc new file mode 100644 index 0000000000..a59a34fd6d --- /dev/null +++ b/tests/simplex/quadrature_lib_01.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test Simplex::PGauss: output its quadrature points and weights. + + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test(const unsigned int n_points) +{ + Simplex::PGauss quad(n_points); + + for (unsigned int q = 0; q < quad.size(); ++q) + { + deallog << quad.point(q) << " "; + deallog << quad.weight(q) << " "; + deallog << std::endl; + } +} + +int +main() +{ + initlog(); + + { + deallog.push("2d-1"); + test<2>(1 /*n_points*/); + deallog.pop(); + } + { + deallog.push("2d-3"); + test<2>(3); + deallog.pop(); + } + { + deallog.push("2d-7"); + test<2>(7); + deallog.pop(); + } + + { + deallog.push("3d-1"); + test<3>(1); + deallog.pop(); + } + { + deallog.push("3d-4"); + test<3>(4); + deallog.pop(); + } + { + deallog.push("3d-10"); + test<3>(10); + deallog.pop(); + } +} diff --git a/tests/simplex/quadrature_lib_01.output b/tests/simplex/quadrature_lib_01.output new file mode 100644 index 0000000000..03258812b2 --- /dev/null +++ b/tests/simplex/quadrature_lib_01.output @@ -0,0 +1,27 @@ + +DEAL:2d-1::0.333333 0.333333 1.00000 +DEAL:2d-3::0.666667 0.166667 0.166667 +DEAL:2d-3::0.166667 0.666667 0.166667 +DEAL:2d-3::0.166667 0.166667 0.166667 +DEAL:2d-7::0.333333 0.333333 0.112500 +DEAL:2d-7::0.797427 0.101287 0.0629696 +DEAL:2d-7::0.101287 0.797427 0.0629696 +DEAL:2d-7::0.101287 0.101287 0.0629696 +DEAL:2d-7::0.0597159 0.470142 0.0661971 +DEAL:2d-7::0.470142 0.0597159 0.0661971 +DEAL:2d-7::0.470142 0.470142 0.0661971 +DEAL:3d-1::0.250000 0.250000 0.250000 0.166667 +DEAL:3d-4::0.138197 0.138197 0.138197 0.0416667 +DEAL:3d-4::0.585410 0.138197 0.138197 0.0416667 +DEAL:3d-4::0.138197 0.585410 0.138197 0.0416667 +DEAL:3d-4::0.138197 0.138197 0.585410 0.0416667 +DEAL:3d-10::0.568431 0.143856 0.143856 0.0362942 +DEAL:3d-10::0.143856 0.143856 0.143856 0.0362942 +DEAL:3d-10::0.143856 0.143856 0.568431 0.0362942 +DEAL:3d-10::0.143856 0.568431 0.143856 0.0362942 +DEAL:3d-10::0.00000 0.500000 0.500000 0.00358166 +DEAL:3d-10::0.500000 0.00000 0.500000 0.00358166 +DEAL:3d-10::0.500000 0.500000 0.00000 0.00358166 +DEAL:3d-10::0.500000 0.00000 0.00000 0.00358166 +DEAL:3d-10::0.00000 0.500000 0.00000 0.00358166 +DEAL:3d-10::0.00000 0.00000 0.500000 0.00358166