/**
* Integration rule for simplex entities.
*
- * Following number of quadrature points are currently supported for 2D and
- * 3D:
+ * Users specify a number `n_points_1D` as an indication of what polynomial
+ * degree to be integrated exactly, similarly to the number of points in a
+ * QGauss quadrature object, even though the present quadrature formula is not
+ * a tensor product. The given value is translated for n_points_1D=1,2,3 to
+ * following number of quadrature points for 2D and 3D:
* - 2D: 1, 3, 7
* - 3D: 1, 4, 10
*
- * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points)`.
+ * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`.
*
* @ingroup simplex
*/
{
public:
/**
- * Constructor taking the number of quadrature points @p n_points.
+ * Constructor taking the number of quadrature points in 1D direction
+ * @p n_points_1D.
*/
- explicit QGauss(const unsigned int n_points);
+ explicit QGauss(const unsigned int n_points_1D);
};
} // namespace Simplex
namespace Simplex
{
template <int dim>
- QGauss<dim>::QGauss(const unsigned int n_points)
+ QGauss<dim>::QGauss(const unsigned int n_points_1D)
: QSimplex<dim>(Quadrature<dim>())
{
// fill quadrature points and quadrature weights
if (dim == 1)
{
- const dealii::QGauss<dim> quad(n_points);
+ const dealii::QGauss<dim> quad(n_points_1D);
this->quadrature_points = quad.get_points();
this->weights = quad.get_weights();
}
else if (dim == 2)
{
- if (n_points == 1)
+ if (n_points_1D == 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)
+ else if (n_points_1D == 2)
{
const double Q23 = 2.0 / 3.0;
const double Q16 = 1.0 / 6.0;
this->weights.emplace_back(Q16);
this->weights.emplace_back(Q16);
}
- else if (n_points == 7)
+ else if (n_points_1D == 3)
{
const double q12 = 0.5;
}
else if (dim == 3)
{
- if (n_points == 1)
+ if (n_points_1D == 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)
+ else if (n_points_1D == 2)
{
const double Q124 = 1.0 / 6.0 / 4.0;
this->weights.emplace_back(Q124);
this->weights.emplace_back(Q124);
}
- else if (n_points == 10)
+ else if (n_points_1D == 3)
{
const double Q16 = 1.0 / 6.0;
{
initlog();
- test(Simplex::FE_P<2>(1), Simplex::QGauss<2>(3));
- test(Simplex::FE_P<2>(2), Simplex::QGauss<2>(7));
- test(Simplex::FE_P<3>(1), Simplex::QGauss<3>(4));
- test(Simplex::FE_P<3>(2), Simplex::QGauss<3>(10));
- test(Simplex::FE_DGP<2>(1), Simplex::QGauss<2>(3));
- test(Simplex::FE_DGP<2>(2), Simplex::QGauss<2>(7));
- test(Simplex::FE_DGP<3>(1), Simplex::QGauss<3>(4));
- test(Simplex::FE_DGP<3>(2), Simplex::QGauss<3>(10));
+ test(Simplex::FE_P<2>(1), Simplex::QGauss<2>(2));
+ test(Simplex::FE_P<2>(2), Simplex::QGauss<2>(3));
+ test(Simplex::FE_P<3>(1), Simplex::QGauss<3>(2));
+ test(Simplex::FE_P<3>(2), Simplex::QGauss<3>(3));
+ test(Simplex::FE_DGP<2>(1), Simplex::QGauss<2>(2));
+ test(Simplex::FE_DGP<2>(2), Simplex::QGauss<2>(3));
+ test(Simplex::FE_DGP<3>(1), Simplex::QGauss<3>(2));
+ test(Simplex::FE_DGP<3>(2), Simplex::QGauss<3>(3));
}
// 3) Select components
Simplex::FE_P<dim> fe(params.degree);
- Simplex::QGauss<dim> quad(dim == 2 ? (params.degree == 1 ? 3 : 7) :
- (params.degree == 1 ? 4 : 10));
+ Simplex::QGauss<dim> quad(params.degree + 1);
- Simplex::QGauss<dim - 1> face_quad(dim == 2 ? (params.degree == 1 ? 2 : 3) :
- (params.degree == 1 ? 3 : 7));
+ Simplex::QGauss<dim - 1> face_quad(params.degree + 1);
Simplex::FE_P<dim> fe_mapping(1);
MappingFE<dim> mapping(fe_mapping);
false,
new Simplex::FE_DGP<dim>(degree),
new MappingFE<dim>(Simplex::FE_P<dim>(1)),
- new Simplex::QGauss<dim>(dim == 2 ? (degree == 1 ? 3 : 7) :
- (degree == 1 ? 4 : 10)),
- new Simplex::QGauss<dim - 1>(dim == 2 ? (degree == 1 ? 2 : 3) :
- (degree == 1 ? 3 : 7)),
+ new Simplex::QGauss<dim>(degree + 1),
+ new Simplex::QGauss<dim - 1>(degree + 1),
initial_refinement,
number_refinement);
}
face_orientation,
false,
false,
- n_points);
- q < n_points;
+ quad_ref.size());
+ q < quad_ref.size();
++q, ++i)
{
deallog << quad.point(i) << " ";
face_orientation,
face_flip,
face_rotation,
- n_points);
- q < n_points;
+ quad_ref.size());
+ q < quad_ref.size();
++q, ++i)
{
deallog << quad.point(i) << " ";
}
{
deallog.push("3d-4");
- test<3>(3);
+ test<3>(2);
deallog.pop();
}
{
deallog.push("3d-10");
- test<3>(7);
+ test<3>(3);
deallog.pop();
}
}
}
{
deallog.push("1d-3");
- test<1>(2);
+ test<1>(3);
deallog.pop();
}
{
}
{
deallog.push("2d-3");
- test<2>(3);
+ test<2>(2);
deallog.pop();
}
{
deallog.push("2d-7");
- test<2>(7);
+ test<2>(3);
deallog.pop();
}
}
{
deallog.push("3d-4");
- test<3>(4);
+ test<3>(2);
deallog.pop();
}
{
deallog.push("3d-10");
- test<3>(10);
+ test<3>(3);
deallog.pop();
}
}
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:1d-3::0.112702 0.277778
+DEAL:1d-3::0.500000 0.444444
+DEAL:1d-3::0.887298 0.277778
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
QGauss<dim - 1> face_quad(degree + 1);
#else
- Simplex::QGauss<dim> quad(dim == 2 ? (degree == 1 ? 3 : 7) :
- (degree == 1 ? 4 : 10));
+ Simplex::QGauss<dim> quad(degree + 1);
- Simplex::QGauss<dim - 1> face_quad(dim == 2 ? (degree == 1 ? 2 : 3) :
- (degree == 1 ? 3 : 7));
+ Simplex::QGauss<dim - 1> face_quad(degree + 1);
#endif
ScratchData<dim> scratch_data(mapping, fe, quad, face_quad);
: triangulation()
, fe(Simplex::FE_P<dim>(degree), dim)
, dof_handler(triangulation)
- , quadrature_formula(Simplex::QGauss<dim>(fe.degree == 1 ? 4 : 10))
+ , quadrature_formula(Simplex::QGauss<dim>(fe.degree + 1))
, mapping(Simplex::FE_P<dim>(1))
, present_time(0.0)
, present_timestep(1.0)