--- /dev/null
+New: The Quadrature class has a new boolean is_tensor_product_flag
+to indicate if the represented quadrature formula is in fact a
+tensor product of (identical) 1D formulas. The corresponding 1D object
+can be queried using Quadrature::get_tensor_basis().
+<br>
+(Daniel Arndt, 2017/07/26)
template <class Archive>
void serialize (Archive &ar, const unsigned int version);
+ /**
+ * This quadrature object is a tensor product in case the tensor product
+ * is formed by the first $size^(1/dim)$ quadrature points.
+ * Return if this is the case.
+ */
+ bool is_tensor_product() const;
+
+ /**
+ * In case the quadrature formula is a tensor product, this function
+ * returns the one-dimensional basis object.
+ * Otherwise, calling this function is not allowed.
+ *
+ * Each time this function is called the constructor for the
+ * one-dimensional Quadrature is called and a copy returned.
+ *
+ * The weights for the one-dimensional formula are rescaled
+ * but they are not modified in case the sum is zero or infinite.
+ */
+ Quadrature<1> get_tensor_basis() const;
+
protected:
/**
* List of quadrature points. To be filled by the constructors of derived
* constructors of derived classes.
*/
std::vector<double> weights;
+
+ /**
+ * This quadrature object is a tensor product in case the tensor product
+ * is formed by the first $size^(1/dim)$ quadrature points
+ */
+ bool is_tensor_product_flag;
};
+template <int dim>
+inline
+bool
+Quadrature<dim>::is_tensor_product () const
+{
+ return is_tensor_product_flag;
+}
+
+
+
template <int dim>
template <class Archive>
inline
Quadrature<0>::Quadrature (const unsigned int n_q)
:
quadrature_points (n_q),
- weights (n_q, 0)
+ weights (n_q, 0),
+ is_tensor_product_flag (false)
{}
Quadrature<dim>::Quadrature (const unsigned int n_q)
:
quadrature_points (n_q, Point<dim>()),
- weights (n_q, 0)
+ weights (n_q, 0),
+ is_tensor_product_flag (dim==1)
{}
const std::vector<double> &weights)
:
quadrature_points(points),
- weights(weights)
+ weights(weights),
+ is_tensor_product_flag (dim==1)
{
Assert (weights.size() == points.size(),
ExcDimensionMismatch(weights.size(), points.size()));
Quadrature<dim>::Quadrature (const std::vector<Point<dim> > &points)
:
quadrature_points(points),
- weights(points.size(), std::atof("Inf"))
+ weights(points.size(), std::atof("Inf")),
+ is_tensor_product_flag (dim==1)
{
Assert(weights.size() == points.size(),
ExcDimensionMismatch(weights.size(), points.size()));
Quadrature<dim>::Quadrature (const Point<dim> &point)
:
quadrature_points(std::vector<Point<dim> > (1, point)),
- weights(std::vector<double> (1, 1.))
+ weights(std::vector<double> (1, 1.)),
+ is_tensor_product_flag (dim==1)
{}
const Quadrature<1> &q2)
:
quadrature_points (q1.size() * q2.size()),
- weights (q1.size() * q2.size())
+ weights (q1.size() * q2.size()),
+ is_tensor_product_flag (true)
{
unsigned int present_index = 0;
for (unsigned int i2=0; i2<q2.size(); ++i2)
const Quadrature<1> &q2)
:
quadrature_points (q2.size()),
- weights (q2.size())
+ weights (q2.size()),
+ is_tensor_product_flag (true)
{
unsigned int present_index = 0;
for (unsigned int i2=0; i2<q2.size(); ++i2)
:
Subscriptor(),
// quadrature_points(1),
- weights(1,1.)
+ weights(1,1.),
+ is_tensor_product_flag (false)
{}
:
Subscriptor(),
quadrature_points (Utilities::fixed_power<dim>(q.size())),
- weights (Utilities::fixed_power<dim>(q.size()))
+ weights (Utilities::fixed_power<dim>(q.size())),
+ is_tensor_product_flag(true)
{
Assert (dim <= 3, ExcNotImplemented());
:
Subscriptor(),
quadrature_points (q.quadrature_points),
- weights (q.weights)
+ weights (q.weights),
+ is_tensor_product_flag (q.is_tensor_product_flag)
{}
{
weights = q.weights;
quadrature_points = q.quadrature_points;
+ is_tensor_product_flag = q.is_tensor_product_flag;
return *this;
}
}
+
+template <int dim>
+Quadrature<1>
+Quadrature<dim>::get_tensor_basis () const
+{
+ Assert (this->is_tensor_product_flag == true,
+ ExcMessage("This function only makes sense if "
+ "this object represents a tensor product!"));
+
+ // Just take the first components of the first quadrature points
+ // and rescale the weights accordingly.
+ const unsigned int n_q_points = this->size();
+ const unsigned int n_q_points_1d
+ = static_cast<unsigned int>(std::round(std::pow(n_q_points, 1./dim)));
+ Assert(Utilities::fixed_power<dim>(n_q_points_1d) == n_q_points,
+ ExcInternalError());
+
+ std::vector<Point<1> > q_points_1d (n_q_points_1d);
+ std::vector<double> weights_1d (n_q_points_1d);
+
+ double sum = 0.;
+ for (unsigned int i=0; i<n_q_points_1d; ++i)
+ {
+ sum += this->weight(i);
+ q_points_1d[i](0) = this->point(i)(0);
+ }
+
+ if (!std::isinf(sum) || sum>1.e-10)
+ for (unsigned int i=0; i<n_q_points_1d; ++i)
+ weights_1d[i] = this->weight(i)/sum;
+ else
+ for (unsigned int i=0; i<n_q_points_1d; ++i)
+ weights_1d[i] = this->weight(i);
+
+ return Quadrature<1> (q_points_1d, weights_1d);
+}
+
+
+
//---------------------------------------------------------------------------
template <int dim>
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the type_trait is_tensor_product for all the quadrature classes
+
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature.h>
+
+template <int dim>
+void
+print_is_tensor_product()
+{
+ Quadrature<dim> q_0;
+ deallog << "Quadrature<" << dim << ">: "
+ << q_0.is_tensor_product() << std::endl;
+
+ QIterated<dim> q_1 (QGauss<1>(1), 1);
+ deallog << "QIterated<" << dim << ">: "
+ << q_1.is_tensor_product() << std::endl;
+
+ QGauss<dim> q_2 (1);
+ deallog << "QGauss<" << dim << ">: "
+ << q_2.is_tensor_product() << std::endl;
+
+ QGaussLobatto<dim> q_3 (2);
+ deallog << "QGaussLobatto<" << dim << ">: "
+ << q_3.is_tensor_product() << std::endl;
+
+ QMidpoint<dim> q_4;
+ deallog << "QMidpoint<" << dim << ">: "
+ << q_4.is_tensor_product() << std::endl;
+
+ QSimpson<dim> q_5;
+ deallog << "QSimpson<" << dim << ">: "
+ << q_5.is_tensor_product() << std::endl;
+
+ QTrapez<dim> q_6;
+ deallog << "QTrapez<" << dim << ">: "
+ << q_6.is_tensor_product() << std::endl;
+
+ QMilne<dim> q_7;
+ deallog << "QMilne<" << dim << ">: "
+ << q_7.is_tensor_product() << std::endl;
+
+ QWeddle<dim> q_8;
+ deallog << "QWeddle<" << dim << ">: "
+ << q_8.is_tensor_product() << std::endl;
+
+ QGaussChebyshev<dim> q_9(3);
+ deallog << "QGaussChebyshev<" << dim << ">: "
+ << q_9.is_tensor_product() << std::endl;
+
+ QGaussRadauChebyshev<dim> q_10(1);
+ deallog << "QGaussRadauChebyshev<" << dim << ">: "
+ << q_10.is_tensor_product() << std::endl;
+
+ QGaussLobattoChebyshev<dim> q_11(2);
+ deallog << "QGaussLobattoChebyshev<" << dim << ">: "
+ << q_11.is_tensor_product() << std::endl;
+
+ QSorted<dim> q_16(q_0);
+ deallog << "QSorted<" << dim << ">: "
+ << q_16.is_tensor_product() << std::endl;
+
+ QTelles<dim> q_17(1, Point<dim>());
+ deallog << "QTelles<" << dim << ">: "
+ << q_17.is_tensor_product() << std::endl;
+}
+
+int main()
+{
+ initlog();
+ deallog << std::boolalpha;
+
+ Quadrature<1> q;
+
+ print_is_tensor_product<1>();
+ QAnisotropic<1> q_1(q);
+ deallog << "QAnisotropic<1>: " << q_1.is_tensor_product() << std::endl;
+ QGaussLog<1> q_13(1);
+ deallog << "QGaussLog<1>: " << q_13.is_tensor_product() << std::endl;
+ QGaussLogR<1> q_14(1);
+ deallog << "QGaussLogR<1>: " << q_14.is_tensor_product() << std::endl;
+
+ print_is_tensor_product<2>();
+ QAnisotropic<2> q_2(q, q);
+ deallog << "QAnisotropic<2>: " << q_2.is_tensor_product() << std::endl;
+ QGaussOneOverR<2> q_15(1, Point<2>());
+ deallog << "QGaussOneOverR<2>: " << q_15.is_tensor_product() << std::endl;
+
+ print_is_tensor_product<3>();
+ QAnisotropic<3> q_3(q, q, q);
+ deallog << "QAnisotropic<3>: " << q_3.is_tensor_product() << std::endl;
+}
+
--- /dev/null
+
+DEAL::Quadrature<1>: true
+DEAL::QIterated<1>: true
+DEAL::QGauss<1>: true
+DEAL::QGaussLobatto<1>: true
+DEAL::QMidpoint<1>: true
+DEAL::QSimpson<1>: true
+DEAL::QTrapez<1>: true
+DEAL::QMilne<1>: true
+DEAL::QWeddle<1>: true
+DEAL::QGaussChebyshev<1>: true
+DEAL::QGaussRadauChebyshev<1>: true
+DEAL::QGaussLobattoChebyshev<1>: true
+DEAL::QSorted<1>: true
+DEAL::QTelles<1>: true
+DEAL::QAnisotropic<1>: true
+DEAL::QGaussLog<1>: true
+DEAL::QGaussLogR<1>: true
+DEAL::Quadrature<2>: false
+DEAL::QIterated<2>: true
+DEAL::QGauss<2>: true
+DEAL::QGaussLobatto<2>: true
+DEAL::QMidpoint<2>: true
+DEAL::QSimpson<2>: true
+DEAL::QTrapez<2>: true
+DEAL::QMilne<2>: true
+DEAL::QWeddle<2>: true
+DEAL::QGaussChebyshev<2>: true
+DEAL::QGaussRadauChebyshev<2>: true
+DEAL::QGaussLobattoChebyshev<2>: true
+DEAL::QSorted<2>: false
+DEAL::QTelles<2>: false
+DEAL::QAnisotropic<2>: false
+DEAL::QGaussOneOverR<2>: false
+DEAL::Quadrature<3>: false
+DEAL::QIterated<3>: true
+DEAL::QGauss<3>: true
+DEAL::QGaussLobatto<3>: true
+DEAL::QMidpoint<3>: true
+DEAL::QSimpson<3>: true
+DEAL::QTrapez<3>: true
+DEAL::QMilne<3>: true
+DEAL::QWeddle<3>: true
+DEAL::QGaussChebyshev<3>: true
+DEAL::QGaussRadauChebyshev<3>: true
+DEAL::QGaussLobattoChebyshev<3>: true
+DEAL::QSorted<3>: false
+DEAL::QTelles<3>: false
+DEAL::QAnisotropic<3>: false
projects / dealii.git / commitdiff