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
+tensor product of 1D formulas. The corresponding 1D objects
can be queried using Quadrature::get_tensor_basis().
<br>
(Daniel Arndt, 2017/07/26)
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.
+ * This function returns true if the quadrature object is a tensor product
+ * of one-dimensional formulas.
*/
bool is_tensor_product() const;
/**
* In case the quadrature formula is a tensor product, this function
- * returns the one-dimensional basis object.
+ * returns the one-dimensional basis objects.
* 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;
+ const std::vector<Quadrature<1>> &get_tensor_basis() const;
protected:
/**
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
+ * Indicates if this object represents quadrature formula that is a tensor
+ * product of one-dimensional formulas.
+ * This flag is set if dim==1 or the constructors taking a Quadrature<1>
+ * (and possibly a Quadrature<dim-1> object) is called.
*/
bool is_tensor_product_flag;
+
+ /**
+ * Stores the one-dimensional tensor basis objects in case this object
+ * can be represented by a tensor product.
+ */
+ std::vector<Quadrature<1> > tensor_basis;
};
DEAL_II_NAMESPACE_OPEN
+namespace
+{
+ template <int dim>
+ void
+ set_tensor_basis_1d(const Quadrature<dim> &,
+ std::vector<Quadrature<1>> &)
+ {}
+
+ template <>
+ void
+ set_tensor_basis_1d(const Quadrature<1> &quadrature,
+ std::vector<Quadrature<1>> &tensor_basis)
+ {
+ tensor_basis = std::vector<Quadrature<1> >(1, quadrature);
+ }
+}
+
+
+
template <>
Quadrature<0>::Quadrature (const unsigned int n_q)
:
quadrature_points (n_q, Point<dim>()),
weights (n_q, 0),
is_tensor_product_flag (dim==1)
-{}
+{
+ set_tensor_basis_1d(*this, this->tensor_basis);
+}
{
Assert (weights.size() == points.size(),
ExcDimensionMismatch(weights.size(), points.size()));
+ set_tensor_basis_1d(*this, this->tensor_basis);
}
{
Assert(weights.size() == points.size(),
ExcDimensionMismatch(weights.size(), points.size()));
+ set_tensor_basis_1d(*this, this->tensor_basis);
}
quadrature_points(std::vector<Point<dim> > (1, point)),
weights(std::vector<double> (1, 1.)),
is_tensor_product_flag (dim==1)
-{}
+{
+ set_tensor_basis_1d(*this, this->tensor_basis);
+}
template <>
:
quadrature_points (q1.size() * q2.size()),
weights (q1.size() * q2.size()),
- is_tensor_product_flag (true)
+ is_tensor_product_flag (q1.is_tensor_product())
{
unsigned int present_index = 0;
for (unsigned int i2=0; i2<q2.size(); ++i2)
Assert ((sum>0.999999) && (sum<1.000001), ExcInternalError());
}
#endif
+
+ if (is_tensor_product_flag)
+ {
+ tensor_basis = q1.get_tensor_basis();
+ tensor_basis.push_back(q2);
+ }
}
:
quadrature_points (q2.size()),
weights (q2.size()),
- is_tensor_product_flag (true)
+ is_tensor_product_flag (true),
+ tensor_basis(std::vector<Quadrature<1> >(1,q2))
{
unsigned int present_index = 0;
for (unsigned int i2=0; i2<q2.size(); ++i2)
weights[k] *= q.weight(i2);
++k;
}
+ for (unsigned int i=0; i<dim; ++i)
+ tensor_basis = std::vector<Quadrature<1> >(dim, q);
}
Subscriptor(),
quadrature_points (q.quadrature_points),
weights (q.weights),
- is_tensor_product_flag (q.is_tensor_product_flag)
+ is_tensor_product_flag (q.is_tensor_product_flag),
+ tensor_basis (q.tensor_basis)
{}
weights = q.weights;
quadrature_points = q.quadrature_points;
is_tensor_product_flag = q.is_tensor_product_flag;
+ tensor_basis = q.tensor_basis;
return *this;
}
template <int dim>
-Quadrature<1>
+const std::vector<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!"));
+ Assert (tensor_basis.size()==dim, ExcInternalError());
- // 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);
+ return tensor_basis;
}
this->weights[k++] = qx.weight(k1);
}
Assert (k==this->size(), ExcInternalError());
+ this->is_tensor_product_flag = true;
+ this->tensor_basis.resize(1);
+ this->tensor_basis[0] = qx;
}
template <int dim>
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
const Quadrature<1> &qy)
- : Quadrature<dim>(qx.size()
- *qy.size())
+ : Quadrature<dim>(qx.size()*qy.size())
{
Assert (dim==2, ExcImpossibleInDim(dim));
unsigned int k=0;
this->weights[k++] = qx.weight(k1) * qy.weight(k2);
}
Assert (k==this->size(), ExcInternalError());
+ this->is_tensor_product_flag = true;
+ this->tensor_basis.resize(2);
+ this->tensor_basis[0] = qx;
+ this->tensor_basis[1] = qy;
}
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
const Quadrature<1> &qy,
const Quadrature<1> &qz)
- : Quadrature<dim>(qx.size()
- *qy.size()
- *qz.size())
+ : Quadrature<dim>(qx.size()*qy.size()*qz.size())
{
Assert (dim==3, ExcImpossibleInDim(dim));
unsigned int k=0;
this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
}
Assert (k==this->size(), ExcInternalError());
+ this->is_tensor_product_flag = true;
+ this->tensor_basis.resize(3);
+ this->tensor_basis[0] = qx;
+ this->tensor_basis[1] = qy;
+ this->tensor_basis[2] = qz;
}
Assert (std::fabs(sum_of_weights-1) < 1e-13,
ExcInternalError());
#endif
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->weights[i-1] = w;
this->weights[n-i] = w;
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(0.5 + 0.5*static_cast<double>(points[i]));
this->weights[i] = 0.5*w[i];
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
{
this->quadrature_points[0] = Point<1>(0.5);
this->weights[0] = 1.0;
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(xpts[i]);
this->weights[i] = wts[i];
};
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(xpts[i]);
this->weights[i] = wts[i];
};
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(xpts[i]);
this->weights[i] = wts[i];
};
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(xpts[i]);
this->weights[i] = wts[i];
};
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
this->weights[i] = revert ? w[n-1-i] : w[i];
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
template <>
"factor out the singularity, which is zero at one point."));
this->weights[i] /= denominator;
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
weights[q] = J * weights[q];
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(p[i]);
this->weights[i] = w[i];
}
-
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(p[i]);
this->weights[i] = w[i];
}
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
this->quadrature_points[i] = Point<1>(p[i]);
this->weights[i] = w[i];
}
-
+ this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
}
--- /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 boolean 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
+check_tensor_product(const std::vector<Quadrature<dim> > &quadratures,
+ const std::vector<std::string> &quadrature_names)
+{
+ Assert(false, ExcNotImplemented());
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<1> > &quadratures,
+ const std::vector<std::string> &quadrature_names)
+{
+ for (unsigned int i=0; i<quadratures.size(); ++i)
+ {
+ const Quadrature<1> &quadrature = quadratures[i];
+ if (quadrature.is_tensor_product())
+ {
+ deallog << "1D " << quadrature_names[i];
+ const auto &q_basis = quadrature.get_tensor_basis();
+ AssertThrow(q_basis.size()==1, ExcInternalError());
+ const auto &q_points = quadrature.get_points();
+ const auto &q_weights = quadrature.get_weights();
+ AssertThrow(q_basis[0].size() == q_points.size(),
+ ExcInternalError());
+ for (unsigned int q=0; q<quadrature.size(); ++q)
+ {
+ std::cout << q_points[q] << " " << q_basis[0].get_points()[q] << std::endl;
+ AssertThrow(std::abs((q_points[q]-q_basis[0].get_points()[q]).norm())<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs(q_weights[q]-q_basis[0].get_weights()[q])<1.e-10,
+ ExcInternalError());
+ }
+ deallog << " OK" << std::endl;
+ }
+ }
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<2> > &quadratures,
+ const std::vector<std::string> &quadrature_names)
+{
+ for (unsigned int i=0; i<quadratures.size(); ++i)
+ {
+ const Quadrature<2> &quadrature = quadratures[i];
+ if (quadrature.is_tensor_product())
+ {
+ deallog << "2D " << quadrature_names[i];
+ const auto &q_basis = quadrature.get_tensor_basis();
+ AssertThrow(q_basis.size()==2, ExcInternalError());
+ const auto &q_points = quadrature.get_points();
+ const auto &q_weights = quadrature.get_weights();
+ AssertThrow(q_basis[0].size()*q_basis[1].size() == q_points.size(),
+ ExcInternalError());
+ unsigned int q=0;
+ for (unsigned int q2=0; q2<q_basis[0].size(); ++q2)
+ for (unsigned int q1=0; q1<q_basis[1].size(); ++q1)
+ {
+ AssertThrow(std::abs(q_points[q][0]-q_basis[0].get_points()[q1][0])<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs(q_points[q][1]-q_basis[1].get_points()[q2][0])<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs((q_weights[q]-q_basis[0].get_weights()[q1]*q_basis[1].get_weights()[q2]))<1.e-10,
+ ExcInternalError());
+ ++q;
+ }
+ deallog << " OK" << std::endl;
+ }
+ }
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<3> > &quadratures,
+ const std::vector<std::string> &quadrature_names)
+{
+ for (unsigned int i=0; i<quadratures.size(); ++i)
+ {
+ const Quadrature<3> &quadrature = quadratures[i];
+ if (quadrature.is_tensor_product())
+ {
+ deallog << "3D " << quadrature_names[i];
+ const auto &q_basis = quadrature.get_tensor_basis();
+ AssertThrow(q_basis.size()==3, ExcInternalError());
+ const auto &q_points = quadrature.get_points();
+ const auto &q_weights = quadrature.get_weights();
+ AssertThrow(q_basis[0].size()*q_basis[1].size()*q_basis[2].size()
+ == q_points.size(), ExcInternalError());
+ unsigned int q=0;
+ for (unsigned int q3=0; q3<q_basis[2].size(); ++q3)
+ for (unsigned int q2=0; q2<q_basis[1].size(); ++q2)
+ for (unsigned int q1=0; q1<q_basis[0].size(); ++q1)
+ {
+ AssertThrow(std::abs(q_points[q][0]-q_basis[0].get_points()[q1][0])<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs(q_points[q][1]-q_basis[1].get_points()[q2][0])<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs(q_points[q][2]-q_basis[2].get_points()[q3][0])<1.e-10,
+ ExcInternalError());
+ AssertThrow(std::abs(q_weights[q]-q_basis[0].get_weights()[q1]*q_basis[1].get_weights()[q2]*q_basis[2].get_weights()[q3])<1.e-10,
+ ExcInternalError());
+ ++q;
+ }
+ deallog << " OK" << std::endl;
+ }
+ }
+}
+
+template <int dim>
+void
+fill_quadrature_vector(std::vector<Quadrature<dim> > &quadratures,
+ std::vector<std::string> &quadrature_names)
+{
+ quadratures.push_back(Quadrature<dim>());
+ quadrature_names.push_back("Quadrature");
+
+ quadratures.push_back(QIterated<dim> (QGauss<1>(2), 2));
+ quadrature_names.push_back("QIterated");
+
+ quadratures.push_back(QGauss<dim> (2));
+ quadrature_names.push_back("QGauss");
+
+ quadratures.push_back(QGaussLobatto<dim> (2));
+ quadrature_names.push_back("QGaussLobatto");
+
+ quadratures.push_back(QMidpoint<dim> ());
+ quadrature_names.push_back("QMidPoint");
+
+ quadratures.push_back(QSimpson<dim> ());
+ quadrature_names.push_back("QSimpson");
+
+ quadratures.push_back(QTrapez<dim> ());
+ quadrature_names.push_back("QTrapez");
+
+ quadratures.push_back(QMilne<dim> ());
+ quadrature_names.push_back("QMilne");
+
+ quadratures.push_back(QWeddle<dim> ());
+ quadrature_names.push_back("QWeddle");
+
+ quadratures.push_back(QGaussChebyshev<dim> (3));
+ quadrature_names.push_back("QGaussChebyshev");
+
+ quadratures.push_back(QGaussRadauChebyshev<dim> (2));
+ quadrature_names.push_back("QGaussRadauChebyshev");
+
+ quadratures.push_back(QGaussLobattoChebyshev<dim> (2));
+ quadrature_names.push_back("QGaussLobattoChebyshev");
+
+ quadratures.push_back(QSorted<dim> (Quadrature<dim>()));
+ quadrature_names.push_back("QSorted");
+
+ quadratures.push_back(QTelles<dim> (1, Point<dim>()));
+ quadrature_names.push_back("QTelles");
+}
+
+int main()
+{
+ initlog();
+ deallog << std::boolalpha;
+
+ Quadrature<1> q;
+
+ std::vector<Quadrature<1> > quadratures_1d;
+ std::vector<std::string> quadrature_names_1d;
+ fill_quadrature_vector(quadratures_1d, quadrature_names_1d);
+ quadratures_1d.push_back(QAnisotropic<1>(q));
+ quadrature_names_1d.push_back("QAnisotropic");
+ quadratures_1d.push_back(QGaussLog<1> (1));
+ quadrature_names_1d.push_back("QGaussLog");
+ quadratures_1d.push_back(QGaussLogR<1> (1));
+ quadrature_names_1d.push_back("QGaussLogR");
+ check_tensor_product(quadratures_1d, quadrature_names_1d);
+
+ std::vector<Quadrature<2> > quadratures_2d;
+ std::vector<std::string> quadrature_names_2d;
+ fill_quadrature_vector(quadratures_2d, quadrature_names_2d);
+ quadratures_2d.push_back(QAnisotropic<2>(q,q));
+ quadrature_names_2d.push_back("QAnisotropic");
+ quadratures_2d.push_back(QGaussOneOverR<2> (1, Point<2>()));
+ quadrature_names_2d.push_back("QGaussOneOverR");
+ check_tensor_product(quadratures_2d, quadrature_names_2d);
+
+ std::vector<Quadrature<3> > quadratures_3d;
+ std::vector<std::string> quadrature_names_3d;
+ fill_quadrature_vector(quadratures_3d, quadrature_names_3d);
+ quadratures_3d.push_back(QAnisotropic<3>(q,q,q));
+ quadrature_names_3d.push_back("QAnisotropic");
+ check_tensor_product(quadratures_3d, quadrature_names_3d);
+}
+
--- /dev/null
+
+DEAL::1D Quadrature OK
+DEAL::1D QIterated OK
+DEAL::1D QGauss OK
+DEAL::1D QGaussLobatto OK
+DEAL::1D QMidPoint OK
+DEAL::1D QSimpson OK
+DEAL::1D QTrapez OK
+DEAL::1D QMilne OK
+DEAL::1D QWeddle OK
+DEAL::1D QGaussChebyshev OK
+DEAL::1D QGaussRadauChebyshev OK
+DEAL::1D QGaussLobattoChebyshev OK
+DEAL::1D QSorted OK
+DEAL::1D QTelles OK
+DEAL::1D QAnisotropic OK
+DEAL::1D QGaussLog OK
+DEAL::1D QGaussLogR OK
+DEAL::2D QIterated OK
+DEAL::2D QGauss OK
+DEAL::2D QGaussLobatto OK
+DEAL::2D QMidPoint OK
+DEAL::2D QSimpson OK
+DEAL::2D QTrapez OK
+DEAL::2D QMilne OK
+DEAL::2D QWeddle OK
+DEAL::2D QGaussChebyshev OK
+DEAL::2D QGaussRadauChebyshev OK
+DEAL::2D QGaussLobattoChebyshev OK
+DEAL::2D QTelles OK
+DEAL::2D QAnisotropic OK
+DEAL::3D QIterated OK
+DEAL::3D QGauss OK
+DEAL::3D QGaussLobatto OK
+DEAL::3D QMidPoint OK
+DEAL::3D QSimpson OK
+DEAL::3D QTrapez OK
+DEAL::3D QMilne OK
+DEAL::3D QWeddle OK
+DEAL::3D QGaussChebyshev OK
+DEAL::3D QGaussRadauChebyshev OK
+DEAL::3D QGaussLobattoChebyshev OK
+DEAL::3D QTelles OK
+DEAL::3D QAnisotropic OK
-// check the type_trait is_tensor_product for all the quadrature classes
+// check the boolean is_tensor_product for all the quadrature classes
#include "../tests.h"
DEAL::QGaussRadauChebyshev<2>: true
DEAL::QGaussLobattoChebyshev<2>: true
DEAL::QSorted<2>: false
-DEAL::QTelles<2>: false
-DEAL::QAnisotropic<2>: false
+DEAL::QTelles<2>: true
+DEAL::QAnisotropic<2>: true
DEAL::QGaussOneOverR<2>: false
DEAL::Quadrature<3>: false
DEAL::QIterated<3>: true
DEAL::QGaussRadauChebyshev<3>: true
DEAL::QGaussLobattoChebyshev<3>: true
DEAL::QSorted<3>: false
-DEAL::QTelles<3>: false
-DEAL::QAnisotropic<3>: false
+DEAL::QTelles<3>: true
+DEAL::QAnisotropic<3>: true
projects / dealii.git / commitdiff