From: Peter Munch Date: Fri, 10 Jul 2020 12:01:31 +0000 (+0200) Subject: Implement Simplex::ScalarPolynomial and Simplex::PGauss also for 1D X-Git-Tag: v9.3.0-rc1~1273^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=56e919077ca7117874c7171af8a0f566bbf9543d;p=dealii.git Implement Simplex::ScalarPolynomial and Simplex::PGauss also for 1D Reorder shape functions --- diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h index 07c3d999e4..02f62ba434 100644 --- a/include/deal.II/simplex/polynomials.h +++ b/include/deal.II/simplex/polynomials.h @@ -40,8 +40,6 @@ namespace Simplex template class ScalarPolynomial : public ScalarPolynomialsBase { - static_assert(dim == 2 || dim == 3, "Dimension not supported!"); - public: /** * Make the dimension available to the outside. diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h index 82c9017d4c..7fc793d4eb 100644 --- a/include/deal.II/simplex/quadrature_lib.h +++ b/include/deal.II/simplex/quadrature_lib.h @@ -28,11 +28,14 @@ namespace Simplex /** * Integration rule for simplex entities. * - * Following number of quadrature points are currently supported: + * Following number of quadrature points are currently supported for 2D and + * 3D: * - 2D: 1, 3, 7 * - 3D: 1, 4, 10 * - * @ingroup simplex + * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points)`. + * + * @ingroup simplex */ template class PGauss : public QSimplex diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc index c1823cfdff..5c33c826dc 100644 --- a/source/simplex/polynomials.cc +++ b/source/simplex/polynomials.cc @@ -25,7 +25,14 @@ namespace Simplex unsigned int compute_n_polynomials(const unsigned int dim, const unsigned int degree) { - if (dim == 2) + if (dim == 1) + { + if (degree == 1) + return 2; + if (degree == 2) + return 3; + } + else if (dim == 2) { if (degree == 1) return 3; @@ -60,22 +67,41 @@ namespace Simplex ScalarPolynomial::compute_value(const unsigned int i, const Point & p) const { - if (dim == 2) + if (dim == 1) { if (this->degree() == 1) { if (i == 0) + return 1.0 - p[0]; + else if (i == 1) return p[0]; + } + else if (this->degree() == 2) + { + if (i == 0) + return 2.0 * p[0] * p[0] - 3.0 * p[0] + 1; else if (i == 1) - return p[1]; + return 2.0 * p[0] * p[0] - p[0]; else if (i == 2) + return -4.0 * p[0] * p[0] + 4.0 * p[0]; + } + } + else if (dim == 2) + { + if (this->degree() == 1) + { + if (i == 0) return 1.0 - p[0] - p[1]; + else if (i == 1) + return p[0]; + else if (i == 2) + return p[1]; } else if (this->degree() == 2) { - const double t1 = p[0]; - const double t2 = p[1]; - const double t3 = 1.0 - p[0] - p[1]; + const double t1 = 1.0 - p[0] - p[1]; + const double t2 = p[0]; + const double t3 = p[1]; if (i == 0) return t1 * (2.0 * t1 - 1.0); @@ -147,24 +173,47 @@ namespace Simplex { Tensor<1, dim> grad; - if (dim == 2) + if (dim == 1) + { + if (this->degree() == 1) + { + if (i == 0) + grad[0] = -1.0; + else if (i == 1) + grad[0] = 1.0; + } + else if (this->degree() == 2) + { + if (i == 0) + grad[0] = 4.0 * p[0] - 3.0; + else if (i == 1) + grad[0] = 4.0 * p[0] - 1.0; + else if (i == 2) + grad[0] = -8.0 * p[0] + 4.0; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + else if (dim == 2) { if (this->degree() == 1) { if (i == 0) { - grad[0] = +1.0; - grad[1] = +0.0; + grad[0] = -1.0; + grad[1] = -1.0; } else if (i == 1) { - grad[0] = +0.0; - grad[1] = +1.0; + grad[0] = +1.0; + grad[1] = +0.0; } else if (i == 2) { - grad[0] = -1.0; - grad[1] = -1.0; + grad[0] = +0.0; + grad[1] = +1.0; } else { @@ -173,32 +222,32 @@ namespace Simplex } else if (this->degree() == 2) { - if (i == 2) + if (i == 0) { grad[0] = -3.0 + 4.0 * (p[0] + p[1]); grad[1] = -3.0 + 4.0 * (p[0] + p[1]); } - else if (i == 0) + else if (i == 1) { grad[0] = 4.0 * p[0] - 1.0; grad[1] = 0.0; } - else if (i == 1) + else if (i == 2) { grad[0] = 0.0; grad[1] = 4.0 * p[1] - 1.0; } - else if (i == 5) + else if (i == 3) { grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]); grad[1] = -4.0 * p[0]; } - else if (i == 3) + else if (i == 4) { grad[0] = 4.0 * p[1]; grad[1] = 4.0 * p[0]; } - else if (i == 4) + else if (i == 5) { grad[0] = -4.0 * p[1]; grad[1] = 4.0 * (1.0 - p[0] - 2.0 * p[1]); @@ -436,6 +485,7 @@ namespace Simplex return std::make_unique>(*this); } + template class ScalarPolynomial<1>; template class ScalarPolynomial<2>; template class ScalarPolynomial<3>; diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc index 94bd3ad9c6..b0ffad033b 100644 --- a/source/simplex/quadrature_lib.cc +++ b/source/simplex/quadrature_lib.cc @@ -33,8 +33,14 @@ namespace Simplex : QSimplex(Quadrature()) { // fill quadrature points and quadrature weights + if (dim == 1) + { + const dealii::QGauss quad(n_points); - if (dim == 2) + this->quadrature_points = quad.get_points(); + this->weights = quad.get_weights(); + } + else if (dim == 2) { if (n_points == 1) { @@ -140,6 +146,7 @@ namespace Simplex } // namespace Simplex +template class Simplex::PGauss<1>; template class Simplex::PGauss<2>; template class Simplex::PGauss<3>; diff --git a/tests/simplex/polynomials_01.cc b/tests/simplex/polynomials_01.cc index 86c276ef9e..eb8c82cfcc 100644 --- a/tests/simplex/polynomials_01.cc +++ b/tests/simplex/polynomials_01.cc @@ -63,6 +63,16 @@ main() { initlog(); + { + deallog.push("1d-1"); + test<1>(1); + deallog.pop(); + } + { + deallog.push("1d-2"); + test<1>(2); + deallog.pop(); + } { deallog.push("2d-1"); test<2>(1); diff --git a/tests/simplex/polynomials_01.output b/tests/simplex/polynomials_01.output index fe28ac73ea..5d8db3cf7e 100644 --- a/tests/simplex/polynomials_01.output +++ b/tests/simplex/polynomials_01.output @@ -1,22 +1,32 @@ -DEAL:2d-1::0.211325 0.211325 0.577350 -DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 -DEAL:2d-1::0.788675 0.211325 2.77556e-17 -DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 -DEAL:2d-1::0.211325 0.788675 0.00000 -DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 -DEAL:2d-2::-0.0872983 -0.0872983 0.425403 0.0508067 0.349193 0.349193 -DEAL:2d-2::-0.549193 0.00000 0.00000 -0.549193 -2.09839 -2.09839 0.450807 0.450807 -0.450807 2.64758 2.64758 -0.450807 -DEAL:2d-2::0.00000 -0.0872983 -0.0872983 0.225403 0.174597 0.774597 -DEAL:2d-2::1.00000 0.00000 0.00000 -0.549193 -0.549193 -0.549193 0.450807 2.00000 -0.450807 1.09839 -0.450807 -2.00000 -DEAL:2d-2::0.687298 -0.0872983 1.38778e-17 0.400000 -6.25620e-18 -4.92550e-17 -DEAL:2d-2::2.54919 0.00000 0.00000 -0.549193 1.00000 1.00000 0.450807 3.54919 -0.450807 -0.450807 -3.54919 -3.54919 -DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.225403 0.774597 0.174597 -DEAL:2d-2::-0.549193 0.00000 0.00000 1.00000 -0.549193 -0.549193 2.00000 0.450807 -2.00000 -0.450807 1.09839 -0.450807 -DEAL:2d-2::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 -DEAL:2d-2::1.00000 0.00000 0.00000 1.00000 1.00000 1.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -2.00000 -DEAL:2d-2::-0.0872983 0.687298 0.00000 0.400000 0.00000 0.00000 -DEAL:2d-2::-0.549193 0.00000 0.00000 2.54919 1.00000 1.00000 3.54919 0.450807 -3.54919 -3.54919 -0.450807 -0.450807 +DEAL:1d-1::0.788675 0.211325 +DEAL:1d-1::-1.00000 1.00000 +DEAL:1d-1::0.211325 0.788675 +DEAL:1d-1::-1.00000 1.00000 +DEAL:1d-2::0.687298 -0.0872983 0.400000 +DEAL:1d-2::-2.54919 -0.549193 3.09839 +DEAL:1d-2::0.00000 0.00000 1.00000 +DEAL:1d-2::-1.00000 1.00000 0.00000 +DEAL:1d-2::-0.0872983 0.687298 0.400000 +DEAL:1d-2::0.549193 2.54919 -3.09839 +DEAL:2d-1::0.577350 0.211325 0.211325 +DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 +DEAL:2d-1::2.77556e-17 0.788675 0.211325 +DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 +DEAL:2d-1::0.00000 0.211325 0.788675 +DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 +DEAL:2d-2::0.425403 -0.0872983 -0.0872983 0.349193 0.0508067 0.349193 +DEAL:2d-2::-2.09839 -2.09839 -0.549193 0.00000 0.00000 -0.549193 2.64758 -0.450807 0.450807 0.450807 -0.450807 2.64758 +DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.774597 0.225403 0.174597 +DEAL:2d-2::-0.549193 -0.549193 1.00000 0.00000 0.00000 -0.549193 -0.450807 -2.00000 0.450807 2.00000 -0.450807 1.09839 +DEAL:2d-2::1.38778e-17 0.687298 -0.0872983 -4.92550e-17 0.400000 -6.25620e-18 +DEAL:2d-2::1.00000 1.00000 2.54919 0.00000 0.00000 -0.549193 -3.54919 -3.54919 0.450807 3.54919 -0.450807 -0.450807 +DEAL:2d-2::-0.0872983 -0.0872983 0.00000 0.174597 0.225403 0.774597 +DEAL:2d-2::-0.549193 -0.549193 -0.549193 0.00000 0.00000 1.00000 1.09839 -0.450807 2.00000 0.450807 -2.00000 -0.450807 +DEAL:2d-2::0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 +DEAL:2d-2::1.00000 1.00000 1.00000 0.00000 0.00000 1.00000 -2.00000 -2.00000 2.00000 2.00000 -2.00000 -2.00000 +DEAL:2d-2::0.00000 -0.0872983 0.687298 0.00000 0.400000 0.00000 +DEAL:2d-2::1.00000 1.00000 -0.549193 0.00000 0.00000 2.54919 -0.450807 -0.450807 3.54919 0.450807 -3.54919 -3.54919 DEAL:3d-1::0.366025 0.211325 0.211325 0.211325 DEAL:3d-1::-1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 DEAL:3d-2::0.214315 -0.0872983 -0.0872983 -0.0872983 0.298387 0.0508067 0.298387 0.298387 0.0508067 0.0508067 diff --git a/tests/simplex/quadrature_lib_01.cc b/tests/simplex/quadrature_lib_01.cc index a59a34fd6d..8e9f8c3002 100644 --- a/tests/simplex/quadrature_lib_01.cc +++ b/tests/simplex/quadrature_lib_01.cc @@ -42,9 +42,24 @@ main() { initlog(); + { + deallog.push("1d-1"); + test<1>(1 /*n_points*/); + deallog.pop(); + } + { + deallog.push("1d-2"); + test<1>(2); + deallog.pop(); + } + { + deallog.push("1d-3"); + test<1>(2); + deallog.pop(); + } { deallog.push("2d-1"); - test<2>(1 /*n_points*/); + test<2>(1); deallog.pop(); } { diff --git a/tests/simplex/quadrature_lib_01.output b/tests/simplex/quadrature_lib_01.output index 03258812b2..24a8420bbd 100644 --- a/tests/simplex/quadrature_lib_01.output +++ b/tests/simplex/quadrature_lib_01.output @@ -1,4 +1,9 @@ +DEAL:1d-1::0.500000 1.00000 +DEAL:1d-2::0.211325 0.500000 +DEAL:1d-2::0.788675 0.500000 +DEAL:1d-3::0.211325 0.500000 +DEAL:1d-3::0.788675 0.500000 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