From d425c7d85266e9e543c439319ca43f15ac65a617 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 30 Mar 2023 14:41:09 -0400 Subject: [PATCH] Fix nodal and Gauss-type quadratures for dim = 0. --- include/deal.II/grid/reference_cell.h | 2 +- source/base/quadrature_lib.cc | 2 + source/grid/reference_cell.cc | 7 +- tests/grid/reference_cell_type_03_q.cc | 10 +++ tests/grid/reference_cell_type_03_q.output | 1 + tests/simplex/nodal_type_quadrature_01.cc | 87 +++++++++++++++++++ tests/simplex/nodal_type_quadrature_01.output | 33 +++++++ 7 files changed, 140 insertions(+), 2 deletions(-) create mode 100644 tests/simplex/nodal_type_quadrature_01.cc create mode 100644 tests/simplex/nodal_type_quadrature_01.output diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 8392edacb1..2a70b4d5a0 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -1226,7 +1226,7 @@ ReferenceCell::vertex(const unsigned int v) const case 0: { if (*this == ReferenceCells::Vertex) - return Point(0); + return Point(); break; } case 1: diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 3e28a1e9a0..9218b04ef0 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -2228,9 +2228,11 @@ template class QGaussSimplex<0>; template class QGaussSimplex<1>; template class QGaussSimplex<2>; template class QGaussSimplex<3>; +template class QGaussWedge<0>; template class QGaussWedge<1>; template class QGaussWedge<2>; template class QGaussWedge<3>; +template class QGaussPyramid<0>; template class QGaussPyramid<1>; template class QGaussPyramid<2>; template class QGaussPyramid<3>; diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index d7e7b0e9b3..2eb2aa2d51 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -848,7 +848,12 @@ operator>>(std::istream &in, ReferenceCell &reference_cell) return in; } - +// explicitly instantiate dimension 0 quadrature in addition to the standard +// dimensions +template Quadrature<0> +ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const; +template const Quadrature<0> & +ReferenceCell::get_nodal_type_quadrature() const; #include "reference_cell.inst" diff --git a/tests/grid/reference_cell_type_03_q.cc b/tests/grid/reference_cell_type_03_q.cc index 755f4cd225..d4975490f1 100644 --- a/tests/grid/reference_cell_type_03_q.cc +++ b/tests/grid/reference_cell_type_03_q.cc @@ -66,6 +66,16 @@ main() { initlog(); + { + deallog.push("0D"); + // It doesn't make sense to integrate in 0D, but make sure that + // get_gauss_type_quadrature() still works + deallog << "0D quadrature size: " + << ReferenceCells::Vertex.get_gauss_type_quadrature<0>(1).size() + << std::endl; + deallog.pop(); + } + { deallog.push("1D"); test<1>(ReferenceCells::Line); diff --git a/tests/grid/reference_cell_type_03_q.output b/tests/grid/reference_cell_type_03_q.output index 42bbf99968..dd920473ba 100644 --- a/tests/grid/reference_cell_type_03_q.output +++ b/tests/grid/reference_cell_type_03_q.output @@ -1,4 +1,5 @@ +DEAL:0D::0D quadrature size: 1 DEAL:1D::ReferenceCell: Line DEAL:1D:: computed barycenter = 0.500000 DEAL:1D:: self-reported barycenter = 0.500000 diff --git a/tests/simplex/nodal_type_quadrature_01.cc b/tests/simplex/nodal_type_quadrature_01.cc new file mode 100644 index 0000000000..d6c09a6243 --- /dev/null +++ b/tests/simplex/nodal_type_quadrature_01.cc @@ -0,0 +1,87 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 - 2023 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 ReferenceCell::get_nodal_type_quadrature() + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test(const ReferenceCell reference_cell) +{ + const Quadrature &quad = + reference_cell.template get_nodal_type_quadrature(); + + for (unsigned int q = 0; q < quad.size(); ++q) + { + // There is nothing to print for dim = 0 in the point + if (dim > 0) + deallog << quad.point(q) << ' '; + deallog << quad.weight(q) << ' '; + deallog << std::endl; + } +} + +int +main() +{ + initlog(); + + { + deallog.push("0d"); + deallog.pop(); + } + { + deallog.push("1d"); + test<1>(ReferenceCells::Line); + deallog.pop(); + } + { + deallog.push("2d-1"); + test<2>(ReferenceCells::Triangle); + deallog.pop(); + } + { + deallog.push("2d-2"); + test<2>(ReferenceCells::Quadrilateral); + deallog.pop(); + } + { + deallog.push("3d-1"); + test<3>(ReferenceCells::Tetrahedron); + deallog.pop(); + } + { + deallog.push("3d-2"); + test<3>(ReferenceCells::Pyramid); + deallog.pop(); + } + { + deallog.push("3d-2"); + test<3>(ReferenceCells::Wedge); + deallog.pop(); + } + { + deallog.push("3d-3"); + test<3>(ReferenceCells::Hexahedron); + deallog.pop(); + } +} diff --git a/tests/simplex/nodal_type_quadrature_01.output b/tests/simplex/nodal_type_quadrature_01.output new file mode 100644 index 0000000000..82bb892f5a --- /dev/null +++ b/tests/simplex/nodal_type_quadrature_01.output @@ -0,0 +1,33 @@ + +DEAL:1d::0.00000 inf +DEAL:1d::1.00000 inf +DEAL:2d-1::0.00000 0.00000 inf +DEAL:2d-1::1.00000 0.00000 inf +DEAL:2d-1::0.00000 1.00000 inf +DEAL:2d-2::0.00000 0.00000 inf +DEAL:2d-2::1.00000 0.00000 inf +DEAL:2d-2::0.00000 1.00000 inf +DEAL:2d-2::1.00000 1.00000 inf +DEAL:3d-1::0.00000 0.00000 0.00000 inf +DEAL:3d-1::1.00000 0.00000 0.00000 inf +DEAL:3d-1::0.00000 1.00000 0.00000 inf +DEAL:3d-1::0.00000 0.00000 1.00000 inf +DEAL:3d-2::-1.00000 -1.00000 0.00000 inf +DEAL:3d-2::1.00000 -1.00000 0.00000 inf +DEAL:3d-2::-1.00000 1.00000 0.00000 inf +DEAL:3d-2::1.00000 1.00000 0.00000 inf +DEAL:3d-2::0.00000 0.00000 1.00000 inf +DEAL:3d-2::0.00000 0.00000 0.00000 inf +DEAL:3d-2::1.00000 0.00000 0.00000 inf +DEAL:3d-2::0.00000 1.00000 0.00000 inf +DEAL:3d-2::0.00000 0.00000 1.00000 inf +DEAL:3d-2::1.00000 0.00000 1.00000 inf +DEAL:3d-2::0.00000 1.00000 1.00000 inf +DEAL:3d-3::0.00000 0.00000 0.00000 inf +DEAL:3d-3::1.00000 0.00000 0.00000 inf +DEAL:3d-3::0.00000 1.00000 0.00000 inf +DEAL:3d-3::1.00000 1.00000 0.00000 inf +DEAL:3d-3::0.00000 0.00000 1.00000 inf +DEAL:3d-3::1.00000 0.00000 1.00000 inf +DEAL:3d-3::0.00000 1.00000 1.00000 inf +DEAL:3d-3::1.00000 1.00000 1.00000 inf -- 2.39.5