#define __fe_H
/*---------------------------- fe.h ---------------------------*/
-#include <lac/dfmatrix.h>
-#include <grid/point.h>
#include <base/exceptions.h>
#include <grid/tria.h>
+#include <lac/dfmatrix.h>
-// forward declarations
-template <int dim> class Boundary;
-template <int dim> class FiniteElement;
-template <int dim> class Quadrature;
-
-
-/**
- Provide a set of flags which tells the #FEValues<>::reinit# function, which
- fields are to be updated for each cell. E.g. if you do not need the
- gradients since you want to assemble the mass matrix, you can switch that
- off. By default, all flags are off, i.e. no reinitialization will be done.
-
- A variable of this type has to be passed to the constructor of the
- #FEValues# object. You can select more than one flag by concatenation
- using the #|# (bitwise #or#) operator.
- */
-enum UpdateFields {
- /**
- * Default: update nothing.
- */
- update_default = 0,
- /**
- * Compute quadrature points in real
- * space (not on unit cell).
- */
- update_q_points = 1,
- /**
- * Transform gradients on unit cell to
- * gradients on real cell.
- */
- update_gradients = 2,
- /**
- * Compute jacobian matrices of the
- * transform between unit and real cell
- * in the evaluation points.
- */
- update_jacobians = 4,
- /**
- * Compute the JxW values (Jacobian
- * determinant at the quadrature point
- * times the weight of this point).
- */
- update_JxW_values = 8,
- /**
- * Compute the points on the real cell
- * on which the ansatz functions are
- * located.
- */
- update_ansatz_points = 16
-};
-
-
-
-
-/**
- Represent a finite element evaluated with a specific quadrature rule.
- 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 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.
-
- The unit cell is defined to be the tensor product of the interval $[0,1]$
- in the present number of dimensions. In part of the literature, the convention
- is used that the unit cell be the tensor product of the interval $[-1,1]$,
- which is to distinguished properly.
-
- 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
- these values are the values and gradients of the shape functions at the
- quadrature points on the real and the unit cell, the location of the
- quadrature points on the real and on the unit cell, the weights of the
- quadrature points, the Jacobian matrices of the mapping from the unit to
- the real cell at the quadrature points and so on.
-
- The Jacobian matrix is defined to be
- $$ J_{ij} = {d\xi_i \over dx_j} $$
- where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are
- the coordinates on the real cell.
- This is the form needed to compute the gradient on the real cell from
- the gradient on the unit cell. If we want to transform the area element
- $dx dy$ from the real to the unit cell, we have to take the determinant of
- the inverse matrix, which is the reciprocal value of the determinant of the
- matrix defined above.
-
- The #FEValues# object keeps track of those fields which really need to
- be computed, since the computation of the gradients of the ansatz functions
- on each real cell can be quite an expensive thing if it is not needed. The
- object knows about which fields are needed by the #UpdateFields# object
- passed through the constructor. In debug mode, the accessor functions, which
- return values from the different fields, check whether the required field
- was initialized, thus avoiding use of unitialized data.
- */
-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
- * specified finite element using the
- * quadrature points of the given
- * quadrature rule.
- *
- * This function actually only fills
- * the fields related to the unit face,
- * the fields related to a real face (like
- * gradients, true quadrature points, etc.)
- * need to be initialized using the
- * #reinit# function.
- */
- FEValues (const FiniteElement<dim> &,
- const Quadrature<dim> &,
- const UpdateFields);
-
- /**
- * Return the value of the #i#th shape
- * function at the #j# quadrature point.
- */
- double shape_value (const unsigned int i,
- const unsigned int j) const;
-
- /**
- * Return a pointer to the matrix holding
- * all values of shape functions at all
- * integration points, on the present cell.
- * For the format of this matrix, see the
- * documentation for the matrix itself.
- */
- const dFMatrix & get_shape_values () const;
-
- /**
- * Return the gradient of the #i#th shape
- * function at the #j# quadrature point.
- * If you want to get the derivative in
- * one of the coordinate directions, use
- * the appropriate function of the #Point#
- * class to extract one component. Since
- * only a reference to the gradient's value
- * is returned, there should be no major
- * performance drawback.
- * The function returns the gradient on the
- * real element, not the reference element.
- */
- const Point<dim> & shape_grad (const unsigned int i,
- const unsigned int j) const;
-
- /**
- * Return a pointer to the matrix holding
- * all gradients of shape functions at all
- * integration points, on the present cell.
- * For the format of this matrix, see the
- * documentation for the matrix itself.
- */
- const vector<vector<Point<dim> > > & get_shape_grads () const;
-
- /**
- * Return the position of the #i#th
- * quadrature point in real space.
- */
- const Point<dim> & quadrature_point (const unsigned int i) const;
-
- /**
- * Return a pointer to the vector of
- * quadrature points.
- */
- const vector<Point<dim> > & get_quadrature_points () const;
-
- /**
- * Return the point in real space where
- * the #i#th ansatz function is located
- * (location is in the sense of where it
- * assumes its nominal properties, e.g. at
- * the vertex of a cell, at the center of
- * a line, etc).
- *
- * This function is needed for the
- * interpolation problem: if we want to
- * transfer a continuous function to a
- * finite element function by interpolation
- * we have to take the continuous
- * function's value at the ansatz function
- * locations.
- */
- const Point<dim> & ansatz_point (const unsigned int i) const;
-
- /**
- * Return a pointer to the vector of points
- * denoting the location of the ansatz
- * functions.
- */
- const vector<Point<dim> > & get_ansatz_points () const;
-
- /**
- * Return the Jacobi determinant times
- * the weight of the #i#th quadrature
- * point.
- */
- double JxW (const unsigned int i) const;
-
- /**
- * Return a pointer to the array holding
- * the JxW values at the different
- * quadrature points.
- */
- const vector<double> & get_JxW_values () const;
-
- /**
- * Reinitialize the gradients, Jacobi
- * determinants, etc for the given cell
- * and the given finite element.
- *
- * This function needs a boundary object
- * passed, since this class needs to know
- * how to handle faces which are located
- * on the boundary of the domain. In that
- * case, faces may be curved and the
- * calculation of quadrature points,
- * gradients and the like may need
- * additional effort, depending on the
- * mapping from the unit to the real cell
- * (linear mappings use straight boundary
- * segments, but higher order elements
- * may use other ways.)
- */
- void reinit (const Triangulation<dim>::cell_iterator &,
- const FiniteElement<dim> &,
- const Boundary<dim> &);
-
- /**
- * Exception
- */
- DeclException2 (ExcInvalidIndex,
- int, int,
- << "The index " << arg1
- << " is out of range, it should be less than " << arg2);
- /**
- * Exception
- */
- DeclException0 (ExcAccessToUninitializedField);
- /**
- * Exception
- */
- DeclException0 (ExcCannotInitializeField);
-
- private:
- /**
- * Store the values of the shape functions
- * at the quadrature points. Rows in this
- * matrix denote the values of a single
- * shape function at the different points,
- * columns are for a single point with the
- * different shape functions.
- */
- dFMatrix shape_values;
-
- /**
- * Store the gradients of the shape
- * functions at the quadrature points.
- * Since unfortunately the full matrix
- * classes of DEAL are not templated,
- * we have to store them in an
- * archetypic style.
- *
- * This field is reset each time
- * #reinit# is called and contains the
- * gradients on the real element, rather
- * than on the reference element.
- */
- vector<vector<Point<dim> > > shape_gradients;
-
- /**
- * Store the gradients of the shape
- * functions at the quadrature points on
- * the unit cell.
- * This field is set up upon construction
- * of the object and contains the gradients
- * on the reference element.
- */
- vector<vector<Point<dim> > > unit_shape_gradients;
-
- /**
- * Store an array of the weights of the
- * quadrature points. This array is
- * set up upon construction.
- */
- vector<double> weights;
-
- /**
- * Store an array of weights times the
- * Jacobi determinant at the quadrature
- * points. This function is reset each time
- * #reinit# is called. The Jacobi determinant
- * is actually the reciprocal value of the
- * Jacobi matrices stored in this class,
- * see the general documentation of this
- * class for more information.
- */
- vector<double> JxW_values;
-
- /**
- * Array of quadrature points. This array
- * is set up upon calling #reinit# and
- * contains the quadrature points on the
- * real element, rather than on the
- * reference element.
- */
- vector<Point<dim> > quadrature_points;
-
- /**
- * Array of quadrature points in the unit
- * cell. This array is set up upon
- * construction and contains the quadrature
- * points on the reference element.
- */
- vector<Point<dim> > unit_quadrature_points;
-
- /**
- * Array of points denoting the off-point
- * of the ansatz functions. In real space
- * (no-one seems to need the off-point
- * on the unit cell, so no function is
- * provided for this).
- */
- vector<Point<dim> > ansatz_points;
-
- /**
- * Store the jacobi matrices at the
- * different quadrature points. This field
- * is set each time #reinit# is called.
- */
- vector<dFMatrix> jacobi_matrices;
-
- /**
- * Store which fields are to be updated by
- * the reinit function.
- */
- UpdateFields update_flags;
-};
-
-
-
-
-/**
- Represent a finite element evaluated with a specific quadrature rule.
- 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 on the unit
- face before doing the quadrature itself. Only the Jacobian matrix of
- the transformation from the unit face to the real face and the integration
- points in real space are calculated each time we move on to a new face.
-
- The unit face is defined to be the tensor product of the interval $[0,1]$
- in the present number of dimensions minus one. In part of the literature,
- the convention is used that the unit cell be the tensor product of the
- interval $[-1,1]$, which is to distinguished properly.
-
- This class is very similar to the #FEValues# class; see there for more
- documentation.
- */
-template <int dim>
-class FEFaceValues {
- public:
-
-
-
- /**
- * Number of quadrature points on
- * the face.
- */
- const unsigned int n_quadrature_points;
-
- /**
- * Total number of shape functions
- * on the cell adjacent to this face.
- * This number is not the same as the
- * number of shape functions of which
- * the center is located on the face.
- */
- const unsigned int total_dofs;
-
- /**
- * Constructor. Fill all arrays with the
- * values of the shape functions of the
- * specified finite element using the
- * quadrature points of the given
- * quadrature rule for the face, which
- * has a dimension one less than the
- * cell.
- *
- * This function actually only fills
- * the fields related to the unit face,
- * the fields related to a real face (like
- * gradients, true quadrature points, etc.)
- * need to be initialized using the
- * #reinit# function.
- */
- FEFaceValues (const FiniteElement<dim> &,
- const Quadrature<dim-1> &,
- const UpdateFields);
-
- /**
- * Return the value of the #i#th shape
- * function at the #j# quadrature point.
- */
- double shape_value (const unsigned int i,
- const unsigned int j) const;
-
- /**
- * Return a pointer to the matrix holding
- * all values of shape functions at all
- * integration points, on the present cell.
- * For the format of this matrix, see the
- * documentation for the matrix itself.
- */
- const dFMatrix & get_shape_values () const;
-
- /**
- * Return the gradient of the #i#th shape
- * function at the #j# quadrature point.
- * If you want to get the derivative in
- * one of the coordinate directions, use
- * the appropriate function of the #Point#
- * class to extract one component. Since
- * only a reference to the gradient's value
- * is returned, there should be no major
- * performance drawback.
- * The function returns the gradient on the
- * real element, not the reference element.
- */
- const Point<dim> & shape_grad (const unsigned int i,
- const unsigned int j) const;
-
- /**
- * Return a pointer to the matrix holding
- * all gradients of shape functions at all
- * integration points, on the present cell.
- * For the format of this matrix, see the
- * documentation for the matrix itself.
- */
- const vector<vector<Point<dim> > > & get_shape_grads () const;
-
- /**
- * Return the position of the #i#th
- * quadrature point in real space.
- *
- * For curved boundary cells, using
- * biquadratic or higher mappings
- * of the unit cell to the real cell,
- * these points may not be on the
- * plane submannifold on which the
- * vertices of the face lie.
- */
- const Point<dim> & quadrature_point (const unsigned int i) const;
-
- /**
- * Return a pointer to the vector of
- * quadrature points.
- */
- const vector<Point<dim> > & get_quadrature_points () const;
-
- /**
- * Return the point in real space where
- * the #i#th ansatz function is located
- * (location is in the sense of where it
- * assumes its nominal properties, e.g. at
- * the vertex of a cell, at the center of
- * a line, etc).
- *
- * This function is needed for the
- * interpolation problem: if we want to
- * transfer a continuous function to a
- * finite element function by interpolation
- * we have to take the continuous
- * function's value at the ansatz function
- * locations.
- */
- const Point<dim> & ansatz_point (const unsigned int i) const;
-
- /**
- * Return a pointer to the vector of points
- * denoting the location of the ansatz
- * functions.
- */
- const vector<Point<dim> > & get_ansatz_points () const;
-
- /**
- * Return the Jacobi determinant times
- * the weight of the #i#th quadrature
- * point.
- */
- double JxW (const unsigned int i) const;
-
- /**
- * Return a pointer to the array holding
- * the JxW values at the different
- * quadrature points.
- */
- const vector<double> & get_JxW_values () const;
-
- /**
- * Reinitialize the gradients, Jacobi
- * determinants, etc for the given cell
- * and the given finite element.
- *
- * The constructor needs a boundary object
- * passed, since this class needs to know
- * how to handle faces which are located
- * on the boundary of the domain. In that
- * case, faces may be curved and the
- * calculation of quadrature points,
- * gradients and the like may need
- * additional effort, depending on the
- * mapping from the unit to the real cell
- * (linear mappings use straight boundary
- * segments, but higher order elements
- * may use other ways.)
- */
- void reinit (const Triangulation<dim>::face_iterator &,
- const FiniteElement<dim> &,
- const Boundary<dim> &);
-
- /**
- * Exception
- */
- DeclException2 (ExcInvalidIndex,
- int, int,
- << "The index " << arg1
- << " is out of range, it should be less than " << arg2);
- /**
- * Exception
- */
- DeclException0 (ExcAccessToUninitializedField);
- /**
- * Exception
- */
- DeclException0 (ExcCannotInitializeField);
-
- private:
- /**
- * Store the values of the shape functions
- * at the quadrature points. Rows in this
- * matrix denote the values of a single
- * shape function at the different points,
- * columns are for a single point with the
- * different shape functions.
- */
- dFMatrix shape_values;
-
- /**
- * Store the gradients of the shape
- * functions at the quadrature points.
- * Since unfortunately the full matrix
- * classes of DEAL are not templated,
- * we have to store them in an
- * archetypic style.
- *
- * This field is reset each time
- * #reinit# is called and contains the
- * gradients on the real element, rather
- * than on the reference element.
- */
- vector<vector<Point<dim> > > shape_gradients;
-
- /**
- * Store the gradients of the shape
- * functions at the quadrature points on
- * the unit cell.
- * This field is set up upon construction
- * of the object and contains the gradients
- * on the reference element.
- */
- vector<vector<Point<dim> > > unit_shape_gradients;
-
- /**
- * Store an array of the weights of the
- * quadrature points. This array is
- * set up upon construction.
- */
- vector<double> weights;
-
- /**
- * Store an array of weights times the
- * Jacobi determinant at the quadrature
- * points. This function is reset each time
- * #reinit# is called. The Jacobi determinant
- * is actually the reciprocal value of the
- * Jacobi matrices stored in this class,
- * see the general documentation of this
- * class for more information.
- */
- vector<double> JxW_values;
-
- /**
- * Array of quadrature points. This array
- * is set up upon calling #reinit# and
- * contains the quadrature points on the
- * real element, rather than on the
- * reference element.
- */
- vector<Point<dim> > quadrature_points;
-
- /**
- * Array of quadrature points in the unit
- * cell. This array is set up upon
- * construction and contains the quadrature
- * points on the reference element.
- */
- vector<Point<dim-1> > unit_quadrature_points;
-
- /**
- * Array of points denoting the off-point
- * of the ansatz functions. In real space
- * (no-one seems to need the off-point
- * on the unit cell, so no function is
- * provided for this).
- */
- vector<Point<dim> > ansatz_points;
-
- /**
- * Store the jacobi matrices at the
- * different quadrature points. This field
- * is set each time #reinit# is called.
- */
- vector<dFMatrix> jacobi_matrices;
-
- /**
- * Store a pointer to the object describing
- * the boundary of the domain.
- */
- const Boundary<dim> &boundary;
-
- /**
- * Store which fields are to be updated by
- * the reinit function.
- */
- UpdateFields update_flags;
-};
-
-
-/*------------------------ Inline functions -----------------------------------*/
-
-
-
-template <int dim>
-inline
-const dFMatrix & FEValues<dim>::get_shape_values () const {
- return shape_values;
-};
-
-
-
-
-template <int dim>
-inline
-const vector<vector<Point<dim> > > &
-FEValues<dim>::get_shape_grads () const {
- Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
- return shape_gradients;
-};
-
-
-
-template <int dim>
-inline
-const vector<Point<dim> > &
-FEValues<dim>::get_quadrature_points () const {
- Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
- return quadrature_points;
-};
-
-
-
-template <int dim>
-inline
-const vector<Point<dim> > &
-FEValues<dim>::get_ansatz_points () const {
- Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
- return ansatz_points;
-};
-
-
-
-template <int dim>
-inline
-const vector<double> &
-FEValues<dim>::get_JxW_values () const {
- Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
- return JxW_values;
-};
-
-
-
/*---------------------------- fe.h ---------------------------*/
/* end of #ifndef __fe_H */
--- /dev/null
+/*---------------------------- fe_values.h ---------------------------*/
+/* $Id$ */
+#ifndef __fe_values_H
+#define __fe_values_H
+/*---------------------------- fe_values.h ---------------------------*/
+
+
+#include <lac/dfmatrix.h>
+#include <grid/point.h>
+#include <base/exceptions.h>
+#include <grid/tria.h>
+
+
+
+// forward declarations
+template <int dim> class Boundary;
+template <int dim> class FiniteElement;
+template <int dim> class Quadrature;
+
+
+/**
+ Provide a set of flags which tells the #FEValues<>::reinit# function, which
+ fields are to be updated for each cell. E.g. if you do not need the
+ gradients since you want to assemble the mass matrix, you can switch that
+ off. By default, all flags are off, i.e. no reinitialization will be done.
+
+ A variable of this type has to be passed to the constructor of the
+ #FEValues# object. You can select more than one flag by concatenation
+ using the #|# (bitwise #or#) operator.
+ */
+enum UpdateFields {
+ /**
+ * Default: update nothing.
+ */
+ update_default = 0,
+ /**
+ * Compute quadrature points in real
+ * space (not on unit cell).
+ */
+ update_q_points = 1,
+ /**
+ * Transform gradients on unit cell to
+ * gradients on real cell.
+ */
+ update_gradients = 2,
+ /**
+ * Compute jacobian matrices of the
+ * transform between unit and real cell
+ * in the evaluation points.
+ */
+ update_jacobians = 4,
+ /**
+ * Compute the JxW values (Jacobian
+ * determinant at the quadrature point
+ * times the weight of this point).
+ */
+ update_JxW_values = 8,
+ /**
+ * Compute the points on the real cell
+ * on which the ansatz functions are
+ * located.
+ */
+ update_ansatz_points = 16
+};
+
+
+
+
+/**
+ Represent a finite element evaluated with a specific quadrature rule.
+ 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 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.
+
+ The unit cell is defined to be the tensor product of the interval $[0,1]$
+ in the present number of dimensions. In part of the literature, the convention
+ is used that the unit cell be the tensor product of the interval $[-1,1]$,
+ which is to distinguished properly.
+
+ 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
+ these values are the values and gradients of the shape functions at the
+ quadrature points on the real and the unit cell, the location of the
+ quadrature points on the real and on the unit cell, the weights of the
+ quadrature points, the Jacobian matrices of the mapping from the unit to
+ the real cell at the quadrature points and so on.
+
+ The Jacobian matrix is defined to be
+ $$ J_{ij} = {d\xi_i \over dx_j} $$
+ where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are
+ the coordinates on the real cell.
+ This is the form needed to compute the gradient on the real cell from
+ the gradient on the unit cell. If we want to transform the area element
+ $dx dy$ from the real to the unit cell, we have to take the determinant of
+ the inverse matrix, which is the reciprocal value of the determinant of the
+ matrix defined above.
+
+ The #FEValues# object keeps track of those fields which really need to
+ be computed, since the computation of the gradients of the ansatz functions
+ on each real cell can be quite an expensive thing if it is not needed. The
+ object knows about which fields are needed by the #UpdateFields# object
+ passed through the constructor. In debug mode, the accessor functions, which
+ return values from the different fields, check whether the required field
+ was initialized, thus avoiding use of unitialized data.
+ */
+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
+ * specified finite element using the
+ * quadrature points of the given
+ * quadrature rule.
+ *
+ * This function actually only fills
+ * the fields related to the unit face,
+ * the fields related to a real face (like
+ * gradients, true quadrature points, etc.)
+ * need to be initialized using the
+ * #reinit# function.
+ */
+ FEValues (const FiniteElement<dim> &,
+ const Quadrature<dim> &,
+ const UpdateFields);
+
+ /**
+ * Return the value of the #i#th shape
+ * function at the #j# quadrature point.
+ */
+ double shape_value (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all values of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const dFMatrix & get_shape_values () const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at the #j# quadrature point.
+ * If you want to get the derivative in
+ * one of the coordinate directions, use
+ * the appropriate function of the #Point#
+ * class to extract one component. Since
+ * only a reference to the gradient's value
+ * is returned, there should be no major
+ * performance drawback.
+ * The function returns the gradient on the
+ * real element, not the reference element.
+ */
+ const Point<dim> & shape_grad (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all gradients of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const vector<vector<Point<dim> > > & get_shape_grads () const;
+
+ /**
+ * Return the position of the #i#th
+ * quadrature point in real space.
+ */
+ const Point<dim> & quadrature_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of
+ * quadrature points.
+ */
+ const vector<Point<dim> > & get_quadrature_points () const;
+
+ /**
+ * Return the point in real space where
+ * the #i#th ansatz function is located
+ * (location is in the sense of where it
+ * assumes its nominal properties, e.g. at
+ * the vertex of a cell, at the center of
+ * a line, etc).
+ *
+ * This function is needed for the
+ * interpolation problem: if we want to
+ * transfer a continuous function to a
+ * finite element function by interpolation
+ * we have to take the continuous
+ * function's value at the ansatz function
+ * locations.
+ */
+ const Point<dim> & ansatz_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of points
+ * denoting the location of the ansatz
+ * functions.
+ */
+ const vector<Point<dim> > & get_ansatz_points () const;
+
+ /**
+ * Return the Jacobi determinant times
+ * the weight of the #i#th quadrature
+ * point.
+ */
+ double JxW (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the array holding
+ * the JxW values at the different
+ * quadrature points.
+ */
+ const vector<double> & get_JxW_values () const;
+
+ /**
+ * Reinitialize the gradients, Jacobi
+ * determinants, etc for the given cell
+ * and the given finite element.
+ *
+ * This function needs a boundary object
+ * passed, since this class needs to know
+ * how to handle faces which are located
+ * on the boundary of the domain. In that
+ * case, faces may be curved and the
+ * calculation of quadrature points,
+ * gradients and the like may need
+ * additional effort, depending on the
+ * mapping from the unit to the real cell
+ * (linear mappings use straight boundary
+ * segments, but higher order elements
+ * may use other ways.)
+ */
+ void reinit (const Triangulation<dim>::cell_iterator &,
+ const FiniteElement<dim> &,
+ const Boundary<dim> &);
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The index " << arg1
+ << " is out of range, it should be less than " << arg2);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcAccessToUninitializedField);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcCannotInitializeField);
+
+ private:
+ /**
+ * Store the values of the shape functions
+ * at the quadrature points. Rows in this
+ * matrix denote the values of a single
+ * shape function at the different points,
+ * columns are for a single point with the
+ * different shape functions.
+ */
+ dFMatrix shape_values;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points.
+ * Since unfortunately the full matrix
+ * classes of DEAL are not templated,
+ * we have to store them in an
+ * archetypic style.
+ *
+ * This field is reset each time
+ * #reinit# is called and contains the
+ * gradients on the real element, rather
+ * than on the reference element.
+ */
+ vector<vector<Point<dim> > > shape_gradients;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points on
+ * the unit cell.
+ * This field is set up upon construction
+ * of the object and contains the gradients
+ * on the reference element.
+ */
+ vector<vector<Point<dim> > > unit_shape_gradients;
+
+ /**
+ * Store an array of the weights of the
+ * quadrature points. This array is
+ * set up upon construction.
+ */
+ vector<double> weights;
+
+ /**
+ * Store an array of weights times the
+ * Jacobi determinant at the quadrature
+ * points. This function is reset each time
+ * #reinit# is called. The Jacobi determinant
+ * is actually the reciprocal value of the
+ * Jacobi matrices stored in this class,
+ * see the general documentation of this
+ * class for more information.
+ */
+ vector<double> JxW_values;
+
+ /**
+ * Array of quadrature points. This array
+ * is set up upon calling #reinit# and
+ * contains the quadrature points on the
+ * real element, rather than on the
+ * reference element.
+ */
+ vector<Point<dim> > quadrature_points;
+
+ /**
+ * Array of quadrature points in the unit
+ * cell. This array is set up upon
+ * construction and contains the quadrature
+ * points on the reference element.
+ */
+ vector<Point<dim> > unit_quadrature_points;
+
+ /**
+ * Array of points denoting the off-point
+ * of the ansatz functions. In real space
+ * (no-one seems to need the off-point
+ * on the unit cell, so no function is
+ * provided for this).
+ */
+ vector<Point<dim> > ansatz_points;
+
+ /**
+ * Store the jacobi matrices at the
+ * different quadrature points. This field
+ * is set each time #reinit# is called.
+ */
+ vector<dFMatrix> jacobi_matrices;
+
+ /**
+ * Store which fields are to be updated by
+ * the reinit function.
+ */
+ UpdateFields update_flags;
+};
+
+
+
+
+/**
+ Represent a finite element evaluated with a specific quadrature rule.
+ 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 on the unit
+ face before doing the quadrature itself. Only the Jacobian matrix of
+ the transformation from the unit face to the real face and the integration
+ points in real space are calculated each time we move on to a new face.
+
+ The unit face is defined to be the tensor product of the interval $[0,1]$
+ in the present number of dimensions minus one. In part of the literature,
+ the convention is used that the unit cell be the tensor product of the
+ interval $[-1,1]$, which is to distinguished properly.
+
+ This class is very similar to the #FEValues# class; see there for more
+ documentation.
+ */
+template <int dim>
+class FEFaceValues {
+ public:
+
+
+
+ /**
+ * Number of quadrature points on
+ * the face.
+ */
+ const unsigned int n_quadrature_points;
+
+ /**
+ * Total number of shape functions
+ * on the cell adjacent to this face.
+ * This number is not the same as the
+ * number of shape functions of which
+ * the center is located on the face.
+ */
+ const unsigned int total_dofs;
+
+ /**
+ * Constructor. Fill all arrays with the
+ * values of the shape functions of the
+ * specified finite element using the
+ * quadrature points of the given
+ * quadrature rule for the face, which
+ * has a dimension one less than the
+ * cell.
+ *
+ * This function actually only fills
+ * the fields related to the unit face,
+ * the fields related to a real face (like
+ * gradients, true quadrature points, etc.)
+ * need to be initialized using the
+ * #reinit# function.
+ */
+ FEFaceValues (const FiniteElement<dim> &,
+ const Quadrature<dim-1> &,
+ const UpdateFields);
+
+ /**
+ * Return the value of the #i#th shape
+ * function at the #j# quadrature point.
+ */
+ double shape_value (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all values of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const dFMatrix & get_shape_values () const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at the #j# quadrature point.
+ * If you want to get the derivative in
+ * one of the coordinate directions, use
+ * the appropriate function of the #Point#
+ * class to extract one component. Since
+ * only a reference to the gradient's value
+ * is returned, there should be no major
+ * performance drawback.
+ * The function returns the gradient on the
+ * real element, not the reference element.
+ */
+ const Point<dim> & shape_grad (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all gradients of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const vector<vector<Point<dim> > > & get_shape_grads () const;
+
+ /**
+ * Return the position of the #i#th
+ * quadrature point in real space.
+ *
+ * For curved boundary cells, using
+ * biquadratic or higher mappings
+ * of the unit cell to the real cell,
+ * these points may not be on the
+ * plane submannifold on which the
+ * vertices of the face lie.
+ */
+ const Point<dim> & quadrature_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of
+ * quadrature points.
+ */
+ const vector<Point<dim> > & get_quadrature_points () const;
+
+ /**
+ * Return the point in real space where
+ * the #i#th ansatz function is located
+ * (location is in the sense of where it
+ * assumes its nominal properties, e.g. at
+ * the vertex of a cell, at the center of
+ * a line, etc).
+ *
+ * This function is needed for the
+ * interpolation problem: if we want to
+ * transfer a continuous function to a
+ * finite element function by interpolation
+ * we have to take the continuous
+ * function's value at the ansatz function
+ * locations.
+ */
+ const Point<dim> & ansatz_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of points
+ * denoting the location of the ansatz
+ * functions.
+ */
+ const vector<Point<dim> > & get_ansatz_points () const;
+
+ /**
+ * Return the Jacobi determinant times
+ * the weight of the #i#th quadrature
+ * point.
+ */
+ double JxW (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the array holding
+ * the JxW values at the different
+ * quadrature points.
+ */
+ const vector<double> & get_JxW_values () const;
+
+ /**
+ * Reinitialize the gradients, Jacobi
+ * determinants, etc for the given cell
+ * and the given finite element.
+ *
+ * The constructor needs a boundary object
+ * passed, since this class needs to know
+ * how to handle faces which are located
+ * on the boundary of the domain. In that
+ * case, faces may be curved and the
+ * calculation of quadrature points,
+ * gradients and the like may need
+ * additional effort, depending on the
+ * mapping from the unit to the real cell
+ * (linear mappings use straight boundary
+ * segments, but higher order elements
+ * may use other ways.)
+ */
+ void reinit (const Triangulation<dim>::face_iterator &,
+ const FiniteElement<dim> &,
+ const Boundary<dim> &);
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The index " << arg1
+ << " is out of range, it should be less than " << arg2);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcAccessToUninitializedField);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcCannotInitializeField);
+
+ private:
+ /**
+ * Store the values of the shape functions
+ * at the quadrature points. Rows in this
+ * matrix denote the values of a single
+ * shape function at the different points,
+ * columns are for a single point with the
+ * different shape functions.
+ */
+ dFMatrix shape_values;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points.
+ * Since unfortunately the full matrix
+ * classes of DEAL are not templated,
+ * we have to store them in an
+ * archetypic style.
+ *
+ * This field is reset each time
+ * #reinit# is called and contains the
+ * gradients on the real element, rather
+ * than on the reference element.
+ */
+ vector<vector<Point<dim> > > shape_gradients;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points on
+ * the unit cell.
+ * This field is set up upon construction
+ * of the object and contains the gradients
+ * on the reference element.
+ */
+ vector<vector<Point<dim> > > unit_shape_gradients;
+
+ /**
+ * Store an array of the weights of the
+ * quadrature points. This array is
+ * set up upon construction.
+ */
+ vector<double> weights;
+
+ /**
+ * Store an array of weights times the
+ * Jacobi determinant at the quadrature
+ * points. This function is reset each time
+ * #reinit# is called. The Jacobi determinant
+ * is actually the reciprocal value of the
+ * Jacobi matrices stored in this class,
+ * see the general documentation of this
+ * class for more information.
+ */
+ vector<double> JxW_values;
+
+ /**
+ * Array of quadrature points. This array
+ * is set up upon calling #reinit# and
+ * contains the quadrature points on the
+ * real element, rather than on the
+ * reference element.
+ */
+ vector<Point<dim> > quadrature_points;
+
+ /**
+ * Array of quadrature points in the unit
+ * cell. This array is set up upon
+ * construction and contains the quadrature
+ * points on the reference element.
+ */
+ vector<Point<dim-1> > unit_quadrature_points;
+
+ /**
+ * Array of points denoting the off-point
+ * of the ansatz functions. In real space
+ * (no-one seems to need the off-point
+ * on the unit cell, so no function is
+ * provided for this).
+ */
+ vector<Point<dim> > ansatz_points;
+
+ /**
+ * Store the jacobi matrices at the
+ * different quadrature points. This field
+ * is set each time #reinit# is called.
+ */
+ vector<dFMatrix> jacobi_matrices;
+
+ /**
+ * Store a pointer to the object describing
+ * the boundary of the domain.
+ */
+ const Boundary<dim> &boundary;
+
+ /**
+ * Store which fields are to be updated by
+ * the reinit function.
+ */
+ UpdateFields update_flags;
+};
+
+
+
+
+
+/*------------------------ Inline functions -----------------------------------*/
+
+
+
+template <int dim>
+inline
+const dFMatrix & FEValues<dim>::get_shape_values () const {
+ return shape_values;
+};
+
+
+
+
+template <int dim>
+inline
+const vector<vector<Point<dim> > > &
+FEValues<dim>::get_shape_grads () const {
+ Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
+ return shape_gradients;
+};
+
+
+
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEValues<dim>::get_quadrature_points () const {
+ Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
+ return quadrature_points;
+};
+
+
+
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEValues<dim>::get_ansatz_points () const {
+ Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
+ return ansatz_points;
+};
+
+
+
+template <int dim>
+inline
+const vector<double> &
+FEValues<dim>::get_JxW_values () const {
+ Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
+ return JxW_values;
+};
+
+
+
+
+
+/*---------------------------- fe_values.h ---------------------------*/
+/* end of #ifndef __fe_values_H */
+#endif
+/*---------------------------- fe_values.h ---------------------------*/
-/*------------------------------- FEValues -------------------------------*/
-
-
-template <int dim>
-FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const UpdateFields update_flags) :
- 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_quadrature_points)),
- unit_shape_gradients(fe.total_dofs,
- 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>()),
- ansatz_points (fe.total_dofs, Point<dim>()),
- jacobi_matrices (quadrature.n_quadrature_points,
- dFMatrix(dim,dim)),
- update_flags (update_flags)
-{
- for (unsigned int i=0; i<fe.total_dofs; ++i)
- 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<n_quadrature_points; ++i)
- {
- weights[i] = quadrature.weight(i);
- unit_quadrature_points[i] = quadrature.quad_point(i);
- };
-};
-
-
-
-template <int dim>
-double FEValues<dim>::shape_value (const unsigned int i,
- const unsigned int j) const {
- Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
- Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
-
- return shape_values(i,j);
-};
-
-
-
-template <int dim>
-const Point<dim> &
-FEValues<dim>::shape_grad (const unsigned int i,
- const unsigned int j) const {
- Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
- Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
- Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
-
- return shape_gradients[i][j];
-};
-
-
-
-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));
- Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
-
- return quadrature_points[i];
-};
-
-
-
-template <int dim>
-const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
- Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
- Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
-
- return ansatz_points[i];
-};
-
-
-
-template <int dim>
-double FEValues<dim>::JxW (const unsigned int i) const {
- Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
- Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
-
- return JxW_values[i];
-};
-
-
-
-template <int dim>
-void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
- const FiniteElement<dim> &fe,
- const Boundary<dim> &boundary) {
- // fill jacobi matrices and real
- // quadrature points
- if ((update_flags | update_jacobians) ||
- (update_flags | update_q_points))
- fe.fill_fe_values (cell,
- unit_quadrature_points,
- jacobi_matrices,
- update_flags | update_jacobians,
- ansatz_points,
- update_flags | update_ansatz_points,
- quadrature_points,
- update_flags | update_q_points,
- boundary);
-
- // compute gradients on real element if
- // requested
- if (update_flags | update_gradients)
- {
- Assert (update_flags | update_jacobians, ExcCannotInitializeField());
-
- for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<n_quadrature_points; ++j)
- for (unsigned int s=0; s<dim; ++s)
- {
- shape_gradients[i][j](s) = 0;
-
- // (grad psi)_s =
- // (grad_{\xi\eta})_b J_{bs}
- // with J_{bs}=(d\xi_b)/(dx_s)
- for (unsigned int b=0; b<dim; ++b)
- shape_gradients[i][j](s)
- +=
- unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
- };
- };
-
-
- // compute Jacobi determinants in
- // quadrature points.
- // refer to the general doc for
- // why we take the inverse of the
- // determinant
- if (update_flags | update_JxW_values)
- {
- Assert (update_flags | update_jacobians, ExcCannotInitializeField());
- for (unsigned int i=0; i<n_quadrature_points; ++i)
- JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
- };
-};
-
-
-
-
-
-
/*------------------------------- FiniteElementBase ----------------------*/
/*------------------------------- Explicit Instantiations -------------*/
-template class FEValues<1>;
-template class FEValues<2>;
-
template class FiniteElementBase<1>;
template class FiniteElementBase<2>;
--- /dev/null
+/* $Id$ */
+
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/quadrature.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_boundary.h>
+
+
+
+/*------------------------------- FEValues -------------------------------*/
+
+
+template <int dim>
+FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const UpdateFields update_flags) :
+ 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_quadrature_points)),
+ unit_shape_gradients(fe.total_dofs,
+ 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>()),
+ ansatz_points (fe.total_dofs, Point<dim>()),
+ jacobi_matrices (quadrature.n_quadrature_points,
+ dFMatrix(dim,dim)),
+ update_flags (update_flags)
+{
+ for (unsigned int i=0; i<fe.total_dofs; ++i)
+ 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<n_quadrature_points; ++i)
+ {
+ weights[i] = quadrature.weight(i);
+ unit_quadrature_points[i] = quadrature.quad_point(i);
+ };
+};
+
+
+
+template <int dim>
+double FEValues<dim>::shape_value (const unsigned int i,
+ const unsigned int j) const {
+ Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+ Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+
+ return shape_values(i,j);
+};
+
+
+
+template <int dim>
+const Point<dim> &
+FEValues<dim>::shape_grad (const unsigned int i,
+ const unsigned int j) const {
+ Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+ Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+ Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
+
+ return shape_gradients[i][j];
+};
+
+
+
+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));
+ Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
+
+ return quadrature_points[i];
+};
+
+
+
+template <int dim>
+const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
+ Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
+ Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
+
+ return ansatz_points[i];
+};
+
+
+
+template <int dim>
+double FEValues<dim>::JxW (const unsigned int i) const {
+ Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+ Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
+
+ return JxW_values[i];
+};
+
+
+
+template <int dim>
+void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
+ const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary) {
+ // fill jacobi matrices and real
+ // quadrature points
+ if ((update_flags | update_jacobians) ||
+ (update_flags | update_q_points))
+ fe.fill_fe_values (cell,
+ unit_quadrature_points,
+ jacobi_matrices,
+ update_flags | update_jacobians,
+ ansatz_points,
+ update_flags | update_ansatz_points,
+ quadrature_points,
+ update_flags | update_q_points,
+ boundary);
+
+ // compute gradients on real element if
+ // requested
+ if (update_flags | update_gradients)
+ {
+ Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+
+ for (unsigned int i=0; i<fe.total_dofs; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ for (unsigned int s=0; s<dim; ++s)
+ {
+ shape_gradients[i][j](s) = 0;
+
+ // (grad psi)_s =
+ // (grad_{\xi\eta})_b J_{bs}
+ // with J_{bs}=(d\xi_b)/(dx_s)
+ for (unsigned int b=0; b<dim; ++b)
+ shape_gradients[i][j](s)
+ +=
+ unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
+ };
+ };
+
+
+ // compute Jacobi determinants in
+ // quadrature points.
+ // refer to the general doc for
+ // why we take the inverse of the
+ // determinant
+ if (update_flags | update_JxW_values)
+ {
+ Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+ };
+};
+
+
+
+
+
+/*------------------------------- Explicit Instantiations -------------*/
+
+template class FEValues<1>;
+template class FEValues<2>;
+