/**
- Base class for quadrature formulae in arbitrary dimensions.
+ Base class for quadrature formulae in arbitrary dimensions. This class
+ stores quadrature points and weights on the unit line [0,1], unit
+ square [0,1]x[0,1], etc. This information is used together with
+ objects of the \Ref{FiniteElement} class to compute the values stored
+ in the \Ref{FEValues} objects.
*/
template <int dim>
class Quadrature {
/**
* Number of quadrature points.
*/
- const unsigned int n_quad_points;
+ const unsigned int n_quadrature_points;
/**
* Constructor.
*/
- Quadrature (const unsigned int n_quad_points);
+ Quadrature (const unsigned int n_quadrature_points);
/**
* Return the #i#th quadrature point.
static const double xpts[] = { 0.288675135, 0.71132486 };
static const double wts[] = { 0.5, 0.5 };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
static const double wts[] = { 0.5*W0, 0.5*W1, 0.5*W1, 0.5*W0,
0.5*W0, 0.5*W1, 0.5*W1, 0.5*W0 };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
static const double xpts[] = { G0, G1, G2, G3 };
static const double wts[] = { W0, W1, W1, W0 };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
static const double xpts[] = { G0, G1, G2, G3, G4, G5, G6, G7 };
static const double wts[] = { W0, W1, W2, W3, W3, W2, W1, W0 };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
static const double xpts[] = { 0.0, 0.5, 1.0 };
static const double wts[] = { 1./6., 2./3., 1./6. };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
static const double xpts[] = { 0.0, 1.0 };
static const double wts[] = { 0.5, 0.5 };
- for (unsigned int i=0; i<n_quad_points; ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points.push_back (Point<1> (xpts[i]));
weights.push_back (wts[i]);
This class is an optimization which avoids evaluating the shape functions
at the quadrature points each time a quadrature takes place. Rather, the
values and gradients (and possibly higher order derivatives in future
- versions of this library) are evaluated once and for all before doing the
- quadrature itself.
+ versions of this library) are evaluated once and for all on the unit
+ cell before doing the quadrature itself. Only the Jacobian matrix of
+ the transformation from the unit cell to the real cell and the integration
+ points in real space are calculated each time we move on to a new cell.
Objects of this class store a multitude of different values needed to
do the assemblage steps on real cells rather than on the unit cell. Among
template <int dim>
class FEValues {
public:
+ /**
+ * Number of quadrature points.
+ */
+ const unsigned int n_quadrature_points;
+
+ /**
+ * Total number of shape functions.
+ */
+ const unsigned int total_dofs;
+
/**
* Constructor. Fill all arrays with the
* values of the shape functions of the
const Point<dim> & shape_grad (const unsigned int i,
const unsigned int j) const;
+ /**
+ * Return the position of the #i#th
+ * quadrature point in real space.
+ */
+ const Point<dim> & quadrature_point (const unsigned int i) const;
+
+ /**
+ * Return the Jacobi determinant times
+ * the weight of the #i#th quadrature
+ * point.
+ */
+ double JxW (const unsigned int i) const;
+
/**
* Reinitialize the gradients, Jacobi
* determinants, etc for the given cell
/**
- Define a finite element type.
+ Define a finite element class.
*/
template <int dim>
class FiniteElement;
template <int dim>
FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
const Quadrature<dim> &quadrature) :
- shape_values(fe.total_dofs, quadrature.n_quad_points),
+ n_quadrature_points(quadrature.n_quadrature_points),
+ total_dofs(fe.total_dofs),
+ shape_values(fe.total_dofs, quadrature.n_quadrature_points),
shape_gradients(fe.total_dofs,
- vector<Point<dim> >(quadrature.n_quad_points)),
+ vector<Point<dim> >(quadrature.n_quadrature_points)),
unit_shape_gradients(fe.total_dofs,
- vector<Point<dim> >(quadrature.n_quad_points)),
- weights(quadrature.n_quad_points, 0),
- JxW_values(quadrature.n_quad_points, 0),
- quadrature_points(quadrature.n_quad_points, Point<dim>()),
- unit_quadrature_points(quadrature.n_quad_points, Point<dim>()),
- jacobi_matrices (quadrature.n_quad_points)
+ vector<Point<dim> >(quadrature.n_quadrature_points)),
+ weights(quadrature.n_quadrature_points, 0),
+ JxW_values(quadrature.n_quadrature_points, 0),
+ quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
+ unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
+ jacobi_matrices (quadrature.n_quadrature_points)
{
for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<quadrature.n_quad_points; ++j)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
{
shape_values(i,j) = fe.shape_value(i, quadrature.quad_point(j));
unit_shape_gradients[i][j]
= fe.shape_grad(i, quadrature.quad_point(j));
};
- for (unsigned int i=0; i<weights.size(); ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
{
weights[i] = quadrature.weight(i);
unit_quadrature_points[i] = quadrature.quad_point(i);
+template <int dim>
+const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const {
+ Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+
+ return quadrature_points[i];
+};
+
+
+
+template <int dim>
+double FEValues<dim>::JxW (const unsigned int i) const {
+ Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+
+ return JxW_values[i];
+};
+
+
+
void FEValues<1>::reinit (const Triangulation<1>::cell_iterator &cell,
const FiniteElement<1> &fe) {
const unsigned int dim=1;
quadrature_points);
// compute gradients on real element
for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<quadrature_points.size(); ++j)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
{
shape_gradients[i][j](s) = 0;
// refer to the general doc for
// why we take the inverse of the
// determinant
- for (unsigned int i=0; i<quadrature_points.size(); ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
};
quadrature_points);
// compute gradients on real element
for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<quadrature_points.size(); ++j)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
{
shape_gradients[i][j](s) = 0;
// refer to the general doc for
// why we take the inverse of the
// determinant
- for (unsigned int i=0; i<quadrature_points.size(); ++i)
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
};
// local mesh width
double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
- unsigned int n_points = unit_points.size();
- for (unsigned int i=0; i<n_points; ++i)
+ for (unsigned int i=0; i<points.size(); ++i)
{
jacobians[i](0,0) = 1./h;
points[i] = cell->vertex(0) + h*unit_points[i];