]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Major reorganisation and support for restriction of finite elements to faces.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 27 Apr 1998 09:29:53 +0000 (09:29 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 27 Apr 1998 09:29:53 +0000 (09:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@197 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/include/fe/fe_update_flags.h [new file with mode: 0644]
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/fe/fe_values.cc

index 187bf90e6302355aae3a70d04bff18d075c2c913..188a130e5cbd83eea2c4f95087276424503cdaa4 100644 (file)
@@ -5,20 +5,43 @@
 /*----------------------------   fe.h     ---------------------------*/
 
 #include <base/exceptions.h>
-#include <grid/tria.h>
+#include <grid/point.h>
+#include <grid/dof.h>
 #include <lac/dfmatrix.h>
 
 
 
+template <int dim> class Boundary;
+template <int dim> struct FiniteElementData;
 
 
 
 /**
-  Base class for finite elements in arbitrary dimensions.
+  Dimension dependent data for finite elements. See the #FiniteElementBase#
+  class for more information.
   */
-template <int dim>
-class FiniteElementBase {
-  public:
+struct FiniteElementData<1> {
+                                    /**
+                                     * Number of degrees of freedom on
+                                     * a vertex.
+                                     */
+    const unsigned int dofs_per_vertex;
+
+                                    /** Number of degrees of freedom
+                                     *  on a line.
+                                     */
+    const unsigned int dofs_per_line;
+
+                                    /**
+                                     * Number of degrees of freedom on a
+                                     * face. This information is
+                                     * redundant to some fields in the
+                                     * derived classes but makes
+                                     * writing dimension independant
+                                     * programs easier.
+                                     */
+    const unsigned int dofs_per_face;
+    
                                     /**
                                      * Total number of degrees of freedom
                                      * on a cell. This information is
@@ -35,58 +58,149 @@ class FiniteElementBase {
                                      * #total_dofs# is therefore a good way to
                                      * check if something went wrong.
                                      */
-    FiniteElementBase () :
+    FiniteElementData () :
+                   dofs_per_vertex(0),
+                   dofs_per_line(0),
+                   dofs_per_face(0),
                    total_dofs(0) {};
-    
+
                                     /**
-                                     * Constructor. You have to set the
-                                     * matrices explicitely after calling
-                                     * this base class' constructor.
+                                     * A more useful version to construct
+                                     * an object of this type.
                                      */
-    FiniteElementBase (const unsigned int total_dofs) :
-                   total_dofs(total_dofs) {};
+    FiniteElementData (const unsigned int dofs_per_vertex,
+                      const unsigned int dofs_per_line) :
+                   dofs_per_vertex(dofs_per_vertex),
+                   dofs_per_line(dofs_per_line),
+                   dofs_per_face(dofs_per_vertex),
+                   total_dofs (2*dofs_per_vertex+dofs_per_line) {};
 
                                     /**
-                                     * Copy constructor.
+                                     * Comparison operator. It is not clear to
+                                     * me why we have to declare and implement
+                                     * this one explicitely.
                                      */
-    FiniteElementBase (const FiniteElementBase &f);
+    bool operator == (const FiniteElementData<1> &) const;
+};
+
+
+
 
+/**
+  Dimension dependent data for finite elements. See the #FiniteElementBase#
+  class for more information.
+  */
+struct FiniteElementData<2> {
                                     /**
-                                     * Destructor. Only declared to have a
-                                     * virtual destructor which the compiler
-                                     * wants to have.
+                                     * Number of degrees of freedom on
+                                     * a vertex.
+                                     */
+    const unsigned int dofs_per_vertex;
+
+                                    /** Number of degrees of freedom
+                                     *  on a line.
+                                     */
+    const unsigned int dofs_per_line;
+
+                                    /** Number of degrees of freedom
+                                     *  on a quad.
+                                     */
+    const unsigned int dofs_per_quad;
+
+                                    /**
+                                     * Number of degrees of freedom on a
+                                     * face. This information is
+                                     * redundant to some fields in the
+                                     * derived classes but makes
+                                     * writing dimension independant
+                                     * programs easier.
                                      */
-    virtual ~FiniteElementBase () {};
+    const unsigned int dofs_per_face;
     
+                                    /**
+                                     * Total number of degrees of freedom
+                                     * on a cell. This information is
+                                     * redundant to some fields in the
+                                     * derived classes but makes
+                                     * writing dimension independant
+                                     * programs easier.
+                                     */
+    const unsigned int total_dofs;
+
+                                    /**
+                                     * Default constructor. Constructs an element
+                                     * which is not so useful. Checking
+                                     * #total_dofs# is therefore a good way to
+                                     * check if something went wrong.
+                                     */
+    FiniteElementData () :
+                   dofs_per_vertex(0),
+                   dofs_per_line(0),
+                   dofs_per_quad(0),
+                   dofs_per_face(0),
+                   total_dofs(0) {};
 
                                     /**
-                                     * Return the value of the #i#th shape
-                                     * function at the point #p#. This function
-                                     * should really be pure, but then we could
-                                     * not make copies of a finite element
-                                     * object even if we did not intend to use
-                                     * this function. Therefore, we omit the
-                                     * #=0# signature and implement this function
-                                     * by throwing an exception.
-                                     * #p# is a point on the reference element.
+                                     * A more useful version to construct
+                                     * an object of this type.
                                      */
-    virtual double shape_value (const unsigned int i,
-                               const Point<dim>& p) const;
+    FiniteElementData (const unsigned int dofs_per_vertex,
+                      const unsigned int dofs_per_line,
+                      const unsigned int dofs_per_quad) :
+                   dofs_per_vertex(dofs_per_vertex),
+                   dofs_per_line(dofs_per_line),
+                   dofs_per_quad(dofs_per_quad),
+                   dofs_per_face(2*dofs_per_vertex+
+                                 dofs_per_line),
+                   total_dofs (4*dofs_per_vertex+
+                               4*dofs_per_quad+
+                               dofs_per_line) {};
 
                                     /**
-                                     * Return the gradient of the #i#th shape
-                                     * function at the point #p#. This function
-                                     * should really be pure, but then we could
-                                     * not make copies of a finite element
-                                     * object even if we did not intend to use
-                                     * this function. Therefore, we omit the
-                                     * #=0# signature and implement this function
-                                     * by throwing an exception.
-                                     * #p# is a point on the reference element,
+                                     * Comparison operator. It is not clear to
+                                     * me why we have to declare and implement
+                                     * this one explicitely.
                                      */
-    virtual Point<dim> shape_grad (const unsigned int i,
-                                  const Point<dim>& p) const;
+    bool operator == (const FiniteElementData<2> &) const;
+};
+
+    
+
+
 
+/**
+  Base class for finite elements in arbitrary dimensions. This class provides
+  several fields which describe a specific finite element and which are filled
+  by derived classes. It more or less only offers the fields and access
+  functions which makes it possible to copy finite elements without knowledge
+  of the actual type (linear, quadratic, etc).
+
+  The implementation of this base class is split into two parts: those fields
+  which are not common to all dimensions (#dofs_per_quad# for example are only
+  useful for #dim>=2#) are put into the #FiniteElementData<dim># class which
+  is explicitely specialized for all used dimensions, while those fields which
+  may be formulated in a dimension-independent way are put into the present
+  class.
+
+  The different matrices are initialized with the correct size, such that in
+  the derived (concrete) finite element classes, their entries must only be
+  filled in; no resizing is needed.
+  */
+template <int dim>
+struct FiniteElementBase : public FiniteElementData<dim> {
+  public:
+                                    /**
+                                     * Construct an object of this type.
+                                     * You have to set the
+                                     * matrices explicitely after calling
+                                     * this base class' constructor. For
+                                     * #dim==1#, #dofs_per_quad# must be
+                                     * zero.
+                                     */
+    FiniteElementBase (const unsigned int dofs_per_vertex,
+                      const unsigned int dofs_per_line,
+                      const unsigned int dofs_per_quad=0);
+    
                                     /**
                                      * Return a readonly reference to the
                                      * matrix which describes the transfer of a
@@ -104,7 +218,7 @@ class FiniteElementBase {
     const dFMatrix & prolongate (const unsigned int child) const;
 
                                     /**
-                                     * Return a readinly reference to the
+                                     * Return a readonly reference to the
                                      * matrix which describes the constraints
                                      * at the interface between a refined and
                                      * an unrefined cell.
@@ -115,68 +229,6 @@ class FiniteElementBase {
                                      */
     const dFMatrix & constraints () const;
     
-                                    /**
-                                     * Compute the Jacobian matrix and the
-                                     * quadrature points as well as the ansatz
-                                     * function locations on the real cell in
-                                     * real space from the given cell
-                                     * and the given quadrature points on the
-                                     * unit cell. The Jacobian matrix is to
-                                     * be computed at every quadrature point.
-                                     * This function has to be in the finite
-                                     * element class, since different finite
-                                     * elements need different transformations
-                                     * of the unit cell to a real cell.
-                                     *
-                                     * The computation of the three fields may
-                                     * share some common code, which is why we
-                                     * put it in one function. However, it may
-                                     * not always be necessary to really
-                                     * compute all fields, so there are
-                                     * bool flags which tell the function which
-                                     * of the fields to actually compute.
-                                     *
-                                     * Refer to the documentation of the
-                                     * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
-                                     *
-                                     * It is provided for the finite element
-                                     * class in one space dimension, but for
-                                     * higher dimensions, it depends on the
-                                     * present fe and needs reimplementation
-                                     * by the user.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
-                                     */
-    virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >               &unit_points,
-                                vector<dFMatrix>    &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &ansatz_points,
-                                const bool           compute_ansatz_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const Boundary<dim> &boundary) const;
-
-                                    /**
-                                     * Return the ansatz points this FE has
-                                     * on a face if a cell would have the
-                                     * given face as a side. This function is
-                                     * needed for higher order elements, if
-                                     * we want to use curved boundary
-                                     * approximations. For that reason, a
-                                     * boundary object has to be passed.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
-                                     */
-    virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
-                                    const Boundary<dim>  &boundary,
-                                    vector<Point<dim> >  &ansatz_points) const;
-    
                                     /**
                                      * Comparison operator. We also check for
                                      * equality of the constraint matrix,
@@ -202,18 +254,20 @@ class FiniteElementBase {
                                     /**
                                      * Exception
                                      */
-    DeclException0 (ExcPureFunctionCalled);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotImplemented);
-                                    /**
-                                     * Exception
-                                     */
     DeclException2 (ExcWrongFieldDimension,
                    int, int,
                    << "The field has not the assumed dimension " << arg2
                    << ", but has " << arg1 << " elements.");
+    DeclException2 (ExcWrongInterfaceMatrixSize,
+                   int, int,
+                   << "The interface matrix has a size of " << arg1
+                   << "x" << arg2
+                   << ", which is not reasonable in the present dimension.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInternalError);
+    
   protected:
                                     /**
                                      * Have #N=2^dim# matrices keeping the
@@ -277,153 +331,36 @@ class FiniteElementBase {
 
 
 
-/**
-  Define a finite element class.
-  */
-template <int dim>
-class FiniteElement;
-
-
 
 /**
-  Finite Element in one dimension.
-
-  {\bf Note on extending the finite element class}
-
-  If you want to extend this class (not by derivation, but by adding new
-  elements), you should be aware that it may be copied at places where you
-  don't expect that (e.g. the #DoFHandler# class keeps a copy). You must
-  thus make sure that the copy operator works correctly, in special if
-  pointers are involved, copying by the default copy constructor supplied
-  by the compiler will result in double deletion and maybe access to data
-  elements which are no more valid.
-
-  Consequence: make sure the copy constructor is correct.
-
-  This class should really be a pure one, with functions having the #=0#
-  signature. However, some instances of this class need to be floating around
-  anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values
-  of #dof_per*# available), so we do not make it pure but rather implement
-  those functions which should in fact be pure to throw an error.
-  */
-class FiniteElement<1> : public FiniteElementBase<1> {
-  public:
-                                    /**
-                                     * Number of degrees of freedom on
-                                     * a vertex.
-                                     */
-    const unsigned int dofs_per_vertex;
-
-                                    /** Number of degrees of freedom
-                                     *  on a line.
-                                     */
-    const unsigned int dofs_per_line;
-
-
-                                    /**
-                                     * Default constructor. The base class
-                                     * produces an invalid element.
-                                     */
-    FiniteElement () :
-                   dofs_per_vertex(0),
-                   dofs_per_line(0) {};
-    
-                                    /**
-                                     * Constructor
-                                     */
-    FiniteElement (const unsigned int dofs_per_vertex,
-                  const unsigned int dofs_per_line) :
-                   FiniteElementBase<1> (2*dofs_per_vertex +
-                                         dofs_per_line),
-                   dofs_per_vertex(dofs_per_vertex),
-                   dofs_per_line  (dofs_per_line) {};
-
-                                    /**
-                                     * Copy constructor
-                                     */
-    FiniteElement (const FiniteElement<1> &fe) :
-                   FiniteElementBase<1> (fe),
-                   dofs_per_vertex(fe.dofs_per_vertex),
-                   dofs_per_line  (fe.dofs_per_line) {};
-
-                                    /**
-                                     * Same pseudo-comparison operator
-                                     * as for the base class.
-                                     */
-    bool operator == (const FiniteElement<1> &f) const;
-
-                                    /**
-                                     * Compute the Jacobian matrix and the
-                                     * quadrature points as well as the ansatz
-                                     * function locations on the real cell in
-                                     * real space from the given cell
-                                     * and the given quadrature points on the
-                                     * unit cell. The Jacobian matrix is to
-                                     * be computed at every quadrature point.
-                                     *
-                                     * Refer to the documentation of the
-                                     * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
-                                     *
-                                     * For one dimensional finite elements,
-                                     * these transformations are usually the
-                                     * same, linear ones, so we provide
-                                     * them in the FE<1> base class. You may,
-                                     * however override this implementation
-                                     * if you would like to use finite elements
-                                     * of higher than first order with
-                                     * non-equidistant integration points, e.g.
-                                     * with an exponential dependence from the
-                                     * distance to the origin. The standard
-                                     * implementation distributes the dofs on
-                                     * the line equidistantly.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
-                                     */
-    virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell,
-                                const vector<Point<1> >               &unit_points,
-                                vector<dFMatrix>  &jacobians,
-                                const bool         compute_jacobians,
-                                vector<Point<1> > &ansatz_points,
-                                const bool         compute_ansatz_points,
-                                vector<Point<1> > &q_points,
-                                const bool         compute_q_points,
-                                const Boundary<1> &boundary) const;
-
-                                    /**
-                                     * Return the ansatz points this FE has
-                                     * on a face if a cell would have the
-                                     * given face as a side. This function is
-                                     * needed for higher order elements, if
-                                     * we want to use curved boundary
-                                     * approximations.
-                                     *
-                                     * Question: is this function useful in 1D?
-                                     * At present it is not implemented.
-                                     */
-    virtual void face_ansatz_points (const Triangulation<1>::face_iterator &face,
-                                    const Boundary<1>  &boundary,
-                                    vector<Point<1> >  &ansatz_points) const;
-};
-
-
-
+  Finite Element in any dimension. This class declares the functionality
+  to fill the fields of the #FiniteElementBase# class. Since this is
+  something that depends on the actual finite element, the functions are
+  declared virtual if it is not possible to provide a reasonable standard
+  implementation.
 
 
-/**
-  Finite Element in two dimensions.
+  {\bf Finite Elements in one dimension}
 
+  Finite elements in one dimension need only set the #restriction# and
+  #prolongation# matrices in #FiniteElementBase<1>#. The constructor of
+  this class in one dimension presets the #interface_constraints# matrix
+  by the unit matrix with dimension one. Changing this behaviour in
+  derived classes is generally not a reasonable idea and you risk getting
+  in terrible trouble.
+  
+  
+  {\bf Finite elements in two dimensions}
+  
   In addition to the fields already present in 1D, a constraint matrix
   is needed in case two quads meet at a common line of which one is refined
   once more than the other one. Then there are constraints referring to the
   hanging nodes on that side of the line which is refined. These constraints
   are represented by a $n\times m$-matrix #line_constraints#, where $n$ is the
-  number of degrees of freedom on the refined side (those dofs on the three
-  vertices plus those on the two lines), and $m$ is that of the unrefined side
+  number of degrees of freedom on the refined side (those dofs on the middle
+  vertex plus those on the two lines), and $m$ is that of the unrefined side
   (those dofs on the two vertices plus those on the line). The matrix is thus
-  a rectangular one, being higher than wide.
+  a rectangular one.
 
   The mapping of the dofs onto the indices of the matrix is as follows:
   let $d_v$ be the number of dofs on a vertex, $d_l$ that on a line, then
@@ -443,70 +380,42 @@ class FiniteElement<1> : public FiniteElementBase<1> {
   degree of freedom to other degrees of freedom which are themselves
   constrained. Only one level of indirection is allowed. It is not known
   at the time of this writing whether this is a constraint itself.
-
-  If you want to extend this class (not by derivation, but by adding new
-  elements), see \Ref{FiniteElement<1>}
-
-  This class should really be a pure one, with functions having the #=0#
-  signature. However, some instances of this class need to be floating around
-  anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values
-  of #dof_per*# available), so we do not make it pure but rather implement
-  those functions which should in fact be pure to throw an error.
   */
-class FiniteElement<2> : public FiniteElementBase<2> {
+template <int dim>
+class FiniteElement : public FiniteElementBase<dim> {
   public:
-                                    /**
-                                     * Number of degrees of freedom on
-                                     * a vertex.
-                                     */
-    const unsigned int dofs_per_vertex;
-
-                                    /** Number of degrees of freedom
-                                     *  on a line.
-                                     */
-    const unsigned int dofs_per_line;
-
-                                    /** Number of degrees of freedom
-                                     *  on a quad.
-                                     */
-    const unsigned int dofs_per_quad;
-
-                                    /**
-                                     * Default constructor. The base class
-                                     * produces an invalid element.
-                                     */
-    FiniteElement () :
-                   dofs_per_vertex(0),
-                   dofs_per_line(0),
-                   dofs_per_quad(0) {};
-
                                     /**
                                      * Constructor
                                      */
     FiniteElement (const unsigned int dofs_per_vertex,
                   const unsigned int dofs_per_line,
-                  const unsigned int dofs_per_quad) :
-                   FiniteElementBase<2> (4*dofs_per_vertex +
-                                         4*dofs_per_line +
-                                         dofs_per_quad),
-                   dofs_per_vertex(dofs_per_vertex),
-                   dofs_per_line  (dofs_per_line),
-                   dofs_per_quad  (dofs_per_quad)  {};
+                  const unsigned int dofs_per_quad=0) :
+                   FiniteElementBase<dim> (dofs_per_vertex,
+                                           dofs_per_line,
+                                           dofs_per_quad) {};
 
                                     /**
-                                     * Copy constructor
+                                     * Destructor. Only declared to have a
+                                     * virtual destructor which the compiler
+                                     * wants to have.
                                      */
-    FiniteElement (const FiniteElement<2> &fe) :
-                   FiniteElementBase<2> (fe),
-                   dofs_per_vertex(fe.dofs_per_vertex),
-                   dofs_per_line  (fe.dofs_per_line),
-                   dofs_per_quad  (fe.dofs_per_quad) {};
+    virtual ~FiniteElement () {};
+    
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at the point #p#.
+                                     * #p# is a point on the reference element.
+                                     */
+    virtual double shape_value (const unsigned int i,
+                               const Point<dim>& p) const = 0;
 
                                     /**
-                                     * Same pseudo-comparison operator
-                                     * as for the base class.
+                                     * Return the gradient of the #i#th shape
+                                     * function at the point #p#.
+                                     * #p# is a point on the reference element,
                                      */
-    bool operator == (const FiniteElement<2> &f) const;
+    virtual Point<dim> shape_grad (const unsigned int i,
+                                  const Point<dim>& p) const = 0;
 
                                     /**
                                      * Compute the Jacobian matrix and the
@@ -516,49 +425,234 @@ class FiniteElement<2> : public FiniteElementBase<2> {
                                      * and the given quadrature points on the
                                      * unit cell. The Jacobian matrix is to
                                      * be computed at every quadrature point.
+                                     * This function has to be in the finite
+                                     * element class, since different finite
+                                     * elements need different transformations
+                                     * of the unit cell to a real cell.
+                                     *
+                                     * The computation of the three fields may
+                                     * share some common code, which is why we
+                                     * put it in one function. However, it may
+                                     * not always be necessary to really
+                                     * compute all fields, so there are
+                                     * bool flags which tell the function which
+                                     * of the fields to actually compute.
                                      *
                                      * Refer to the documentation of the
                                      * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
+                                     * of the Jacobi matrix and of the various
+                                     * structures to be filled.
                                      *
-                                     * For two dimensional finite elements,
-                                     * these transformations are usually
-                                     * dependent on the actual finite element,
-                                     * which is expressed by the names
-                                     * sub- and isoparametric elements. This
-                                     * function is therefore not implemented
-                                     * by the FE<2> base class, but is made
-                                     * pure virtual.
+                                     * It is provided for the finite element
+                                     * class in one space dimension, but for
+                                     * higher dimensions, it depends on the
+                                     * present fe and needs reimplementation
+                                     * by the user. This is due to the fact
+                                     * that the user may want to use
+                                     * iso- or subparametric mappings of the
+                                     * unit cell to the real cell, which
+                                     * makes things much more complicated.
                                      *
                                      * The function assumes that the fields
                                      * already have the right number of
                                      * elements.
+                                     *
+                                     * This function is more or less an
+                                     * interface to the #FEValues# class and
+                                     * should not be used by users unless
+                                     * absolutely needed.
                                      */
-    virtual void fill_fe_values (const Triangulation<2>::cell_iterator &cell,
-                                const vector<Point<2> >               &unit_points,
-                                vector<dFMatrix>  &jacobians,
-                                const bool         compute_jacobians,
-                                vector<Point<2> > &ansatz_points,
-                                const bool         compute_ansatz_points,
-                                vector<Point<2> > &q_points,
-                                const bool         compute_q_points,
-                                const Boundary<2> &boundary) const;
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
+                                vector<dFMatrix>    &jacobians,
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points,
+                                const Boundary<dim> &boundary) const;
 
                                     /**
-                                     * Return the ansatz points this FE has
-                                     * on a face if a cell would have the
-                                     * given face as a side. This function is
-                                     * needed for higher order elements, if
-                                     * we want to use curved boundary
-                                     * approximations.
+                                     * Do the same thing that the other
+                                     * #fill_fe_values# function does,
+                                     * exception that a face rather than
+                                     * a cell is considered. The #face_no#
+                                     * parameter denotes the number of the
+                                     * face to the given cell to be
+                                     * considered.
+                                     *
+                                     * The unit points for the quadrature
+                                     * formula are given on the unit face
+                                     * which is a mannifold of dimension
+                                     * one less than the dimension of the
+                                     * cell. The #global_unit_points# 
+                                     * denote the position of the unit points
+                                     * on the selected face on the unit cell.
+                                     * This additional information is passed
+                                     * since the #FEFaceValues# class can
+                                     * compute them once and for all,
+                                     * eliminating the need to recompute it
+                                     * each time #FEFaceValues::reinit# is
+                                     * called.
+                                     *
+                                     * The jacobian matrix is evaluated at
+                                     * each of the quadrature points on the
+                                     * given face. The matrix is the
+                                     * transformation matrix of the unit cell
+                                     * to the real cell, not from the unit
+                                     * face to the real face. This is the
+                                     * necessary matrix to compute the real
+                                     * gradients.
+                                     *
+                                     * Conversely, the Jacobi determinants
+                                     * are the determinants of the
+                                     * transformation from the unit face to
+                                     * the real face. This information is
+                                     * needed to actually perform integrations
+                                     * along faces. Note that we here return
+                                     * the inverse of the determinant of the
+                                     * jacobi matrices as explained in the
+                                     * documentation of the #FEValues# class.
+                                     * 
+                                     * The ansatz points are the
+                                     * off-points of those ansatz functions
+                                     * located on the given face; this
+                                     * information is taken over from the
+                                     * #get_face_ansatz_points# function.
+                                     *
+                                     * The order of ansatz functions is the
+                                     * same as if it were a cell of dimension
+                                     * one less than the present. E.g. in
+                                     * two dimensions, the order is first
+                                     * the vertex functions (using the
+                                     * direction of the face induced by the
+                                     * given cell) then the interior functions.
+                                     * The same applies for the quadrature
+                                     * points which also use the standard
+                                     * direction of faces as laid down by
+                                     * the #Triangulation# class.
+                                     *
+                                     * There is a standard implementation for
+                                     * dimensions greater than one. It
+                                     * uses the #fill_fe_values()#
+                                     * function to retrieve the wanted
+                                     * information. Since this operation acts
+                                     * only on unit faces and cells it does
+                                     * not depend on a specific finite element
+                                     * transformation and is thus applicable
+                                     * for all finite elements and uses tha
+                                     * same mapping from the unit to the real
+                                     * cell as used for the other operations
+                                     * performed by the specific finite element
+                                     * class.
+                                     *
+                                     * Two fields remain to be finite element
+                                     * specific in this standard implementation:
+                                     * The jacobi determinants of the
+                                     * transformation from the unit face to the
+                                     * real face and the ansatz points. For
+                                     * these two fields, there exist pure
+                                     * virtual functions, #get_face_jacobians#
+                                     * and the #get_face_ansatz_points#
+                                     * function.
+                                     *
+                                     * Though there is a standard
+                                     * implementation, there
+                                     * may be room for optimizations which is
+                                     * why this function is made virtual.
+                                     *
+                                     * Since any implementation for one
+                                     * dimension would be senseless, all
+                                     * derived classes should throw an error
+                                     * when called with #dim==1#.
                                      *
                                      * The function assumes that the fields
                                      * already have the right number of
                                      * elements.
+                                     *
+                                     * This function is more or less an
+                                     * interface to the #FEValues# class and
+                                     * should not be used by users unless
+                                     * absolutely needed.
+                                     */
+    virtual void fill_fe_face_values (const DoFHandler<dim>::cell_iterator &cell,
+                                     const unsigned int           face_no,
+                                     const vector<Point<dim-1> > &unit_points,
+                                     const vector<Point<dim> >   &global_unit_points,
+                                     vector<dFMatrix>    &jacobians,
+                                     const bool           compute_jacobians,
+                                     vector<Point<dim> > &ansatz_points,
+                                     const bool           compute_ansatz_points,
+                                     vector<Point<dim> > &q_points,
+                                     const bool           compute_q_points,
+                                     vector<double>      &face_jacobi_determinants,
+                                     const bool           compute_face_jacobians,
+                                     const Boundary<dim> &boundary) const;
+    
+                                    /**
+                                     * This function produces a subset of
+                                     * the information provided by the
+                                     * #fill_fe_face_values()# function,
+                                     * namely the ansatz function off-points
+                                     * of those ansatz functions located on
+                                     * the face. However, you should not try
+                                     * to implement this function using the
+                                     * abovementioned function, since usually
+                                     * that function uses this function to
+                                     * compute information.
+                                     *
+                                     * This function is excluded from the
+                                     * abovementioned one, since no information
+                                     * about the neighboring cell is needed,
+                                     * such that loops over faces alone are
+                                     * possible when using this function.
+                                     * This is useful for example if we want
+                                     * to interpolate boundary values to the
+                                     * finite element functions. If integration
+                                     * along faces is needed, we still need
+                                     * the #fill_fe_values# function.
+                                     *
+                                     * The function assumes that the
+                                     * #ansatz_points# array already has the
+                                     * right size. The order of points in
+                                     * the array matches that returned by
+                                     * the #face->get_dof_indices# function.
+                                     *
+                                     * Since any implementation for one
+                                     * dimension would be senseless, all
+                                     * derived classes should throw an error
+                                     * when called with #dim==1#.
                                      */
-    virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face,
-                                    const Boundary<2>  &boundary,
-                                    vector<Point<2> >  &ansatz_points) const;
+    virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+                                        const Boundary<dim> &boundary,
+                                        vector<Point<dim> > &ansatz_points) const =0;
+
+                                    /**
+                                     * This is the second separated function
+                                     * described in the documentation of the
+                                     * #fill_fe_face_values# function. It
+                                     * returns the determinants of the
+                                     * transformation from the unit face to the
+                                     * real face at the
+                                     *
+                                     * Since any implementation for one
+                                     * dimension would be senseless, all
+                                     * derived classes should throw an error
+                                     * when called with #dim==1#.
+                                     */
+    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<double>      &face_jacobi_determinants) const =0;
+    
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcPureFunctionCalled);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
 };
 
 
index 792afc5073e496c3a27d52dc443b6ba3e5cae835..ab435d964e1bb1baca1a4b528dfa486162cfb735 100644 (file)
@@ -45,36 +45,20 @@ class FELinear : public FiniteElement<dim> {
                                  const Point<dim>& p) const;
 
                                     /**
-                                     * Compute the Jacobian matrix and the
-                                     * quadrature points as well as the ansatz
-                                     * function locations on the real cell in
-                                     * real space from the given cell
-                                     * and the given quadrature points on the
-                                     * unit cell. The Jacobian matrix is to
-                                     * be computed at every quadrature point.
-                                     *
-                                     * Refer to the documentation of the
-                                     * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
+                                     * Refer to the base class for detailed
+                                     * information on this function.
                                      *
                                      * For one dimensional elements, this
                                      * function simply passes through to
                                      * the one implemented in the base class.
-                                     * For two dimensional finite elements,
-                                     * these transformations are usually
-                                     * dependent on the actual finite element,
-                                     * which is expressed by the names
-                                     * sub- and isoparametric elements. This
-                                     * function is therefore not implemented
-                                     * by the FE<2> base class, but is made
-                                     * pure virtual.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
+                                     * For higher dimensional finite elements
+                                     * we use linear mappings and therefore
+                                     * the boundary object is ignored since
+                                     * the boundary is approximated using
+                                     * piecewise straight boundary segments.
                                      */
-    virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >               &unit_points,
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &ansatz_points,
@@ -84,20 +68,21 @@ class FELinear : public FiniteElement<dim> {
                                 const Boundary<dim> &boundary) const;
 
                                     /**
-                                     * Return the ansatz points this FE has
-                                     * on a face if a cell would have the
-                                     * given face as a side. Since we have no
-                                     * degrees of freedom on the faces for
-                                     * the linear ansatz, the ansatz points are
-                                     * simply the vertices of the face.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
+                                     * Refer to the base class for detailed
+                                     * information on this function.
                                      */
-    virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
-                                    const Boundary<dim>  &boundary,
-                                    vector<Point<dim> >  &ansatz_points) const;
+    virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+                                        const Boundary<dim> &boundary,
+                                        vector<Point<dim> > &ansatz_points) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<double>      &face_jacobi_determinants) const;
 };
 
 
@@ -131,36 +116,15 @@ class FEQuadratic : public FiniteElement<dim> {
                                  const Point<dim>& p) const;
 
                                     /**
-                                     * Compute the Jacobian matrix and the
-                                     * quadrature points as well as the ansatz
-                                     * function locations on the real cell in
-                                     * real space from the given cell
-                                     * and the given quadrature points on the
-                                     * unit cell. The Jacobian matrix is to
-                                     * be computed at every quadrature point.
-                                     *
-                                     * Refer to the documentation of the
-                                     * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
+                                     * Refer to the base class for detailed
+                                     * information on this function.
                                      *
                                      * For one dimensional elements, this
                                      * function simply passes through to
                                      * the one implemented in the base class.
-                                     * For two dimensional finite elements,
-                                     * these transformations are usually
-                                     * dependent on the actual finite element,
-                                     * which is expressed by the names
-                                     * sub- and isoparametric elements. This
-                                     * function is therefore not implemented
-                                     * by the FE<2> base class, but is made
-                                     * pure virtual.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
                                      */
-    virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >               &unit_points,
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &ansatz_points,
@@ -168,6 +132,23 @@ class FEQuadratic : public FiniteElement<dim> {
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
                                 const Boundary<dim> &boundary) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+                                        const Boundary<dim> &boundary,
+                                        vector<Point<dim> > &ansatz_points) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<double>      &face_jacobi_determinants) const;
 };
 
 
@@ -201,36 +182,15 @@ class FECubic : public FiniteElement<dim> {
                                  const Point<dim>& p) const;
 
                                     /**
-                                     * Compute the Jacobian matrix and the
-                                     * quadrature points as well as the ansatz
-                                     * function locations on the real cell in
-                                     * real space from the given cell
-                                     * and the given quadrature points on the
-                                     * unit cell. The Jacobian matrix is to
-                                     * be computed at every quadrature point.
-                                     *
-                                     * Refer to the documentation of the
-                                     * \Ref{FEValues} class for a definition
-                                     * of the Jacobi matrix.
+                                     * Refer to the base class for detailed
+                                     * information on this function.
                                      *
                                      * For one dimensional elements, this
                                      * function simply passes through to
                                      * the one implemented in the base class.
-                                     * For two dimensional finite elements,
-                                     * these transformations are usually
-                                     * dependent on the actual finite element,
-                                     * which is expressed by the names
-                                     * sub- and isoparametric elements. This
-                                     * function is therefore not implemented
-                                     * by the FE<2> base class, but is made
-                                     * pure virtual.
-                                     *
-                                     * The function assumes that the fields
-                                     * already have the right number of
-                                     * elements.
                                      */
-    virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >               &unit_points,
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &ansatz_points,
@@ -238,6 +198,23 @@ class FECubic : public FiniteElement<dim> {
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
                                 const Boundary<dim> &boundary) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+                                        const Boundary<dim> &boundary,
+                                        vector<Point<dim> > &ansatz_points) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<double>      &face_jacobi_determinants) const;
 };
 
 
@@ -246,3 +223,5 @@ class FECubic : public FiniteElement<dim> {
 /* end of #ifndef __fe_lib_H */
 #endif
 /*----------------------------   fe_lib.h     ---------------------------*/
+
+
diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h
new file mode 100644 (file)
index 0000000..c6562c8
--- /dev/null
@@ -0,0 +1,59 @@
+/*----------------------------   fe_update_flags.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __fe_update_flags_H
+#define __fe_update_flags_H
+/*----------------------------   fe_update_flags.h     ---------------------------*/
+
+
+/**
+  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 UpdateFlags {
+                                      /**
+                                       * 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
+};
+
+
+
+
+/*----------------------------   fe_update_flags.h     ---------------------------*/
+/* end of #ifndef __fe_update_flags_H */
+#endif
+/*----------------------------   fe_update_flags.h     ---------------------------*/
index caa8d89214feb0dd48175d1be46fd8ab5d702593..e4946a8299b247f38a4f93458c145f3dfd22cf14 100644 (file)
@@ -9,7 +9,7 @@
 #include <grid/point.h>
 #include <base/exceptions.h>
 #include <grid/tria.h>
-
+#include <fe/fe_update_flags.h>
 
 
 // forward declarations
@@ -18,51 +18,6 @@ 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
-};
-
 
 
 
@@ -102,7 +57,7 @@ enum UpdateFields {
   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
+  object knows about which fields are needed by the #UpdateFlags# 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.
@@ -139,7 +94,7 @@ class FEValues {
                                      */
     FEValues (const FiniteElement<dim> &,
              const Quadrature<dim> &,
-             const UpdateFields);
+             const UpdateFlags);
 
                                     /**
                                      * Return the value of the #i#th shape
@@ -251,7 +206,7 @@ class FEValues {
                                      * segments, but higher order elements
                                      * may use other ways.)
                                      */
-    void reinit (const Triangulation<dim>::cell_iterator &,
+    void reinit (const DoFHandler<dim>::cell_iterator &,
                 const FiniteElement<dim> &,
                 const Boundary<dim> &);
 
@@ -363,7 +318,7 @@ class FEValues {
                                      * Store which fields are to be updated by
                                      * the reinit function.
                                      */
-    UpdateFields         update_flags;
+    UpdateFlags         update_flags;
 };
 
 
@@ -385,7 +340,52 @@ class FEValues {
   interval $[-1,1]$, which is to distinguished properly.
 
   This class is very similar to the #FEValues# class; see there for more
-  documentation.
+  documentation. It is, however, a bit more involved: since we want to
+  compute the restriction of finite element functions (here: the basis
+  functions, but a finite element function is obtained by multiplication
+  with the nodal values and summation) to the face of a cell and since
+  finite element functions and especially their gradients need not be
+  continuous at faces, we can not compute the wanted information from
+  the face and a finite element class on the unit cell alone, but we
+  need the real cell as well. In addition, we need to know what number
+  the face is in the set of faces of the cell we want to restrict.
+  Finally, since we may want to use higher order elements with unit cell
+  to real cell mappings of higher than first order, thus applying curved
+  boundaries, we need to know an object describing the boundary of the
+  domain.
+
+
+  {\bf Technical issues}
+
+  Just like in the #FEValues# class, function values and gradients on the unit
+  cell are evaluated at the quadrature points only once, in the constructor.
+  Being a tensor of rank zero, the function values remain the same when we
+  want them at the quadrature points on the real cell, while we get the
+  gradients (a tensor of rank one) by multiplication with the Jacobi matrix
+  of the transformation, which we need to compute for each cell and each
+  quadrature point.
+
+  However, while in the #FEValues# class the quadrature points are always the
+  same, here we deal with more than one face. We therefore store the values
+  and gradients of the ansatz functions on the unit cell in an array with as
+  many elements as there are faces on a cell. The same applies for the
+  quadrature points on the faces: for each face we store the position on the
+  cell. This way we still need to evaluate unit gradients and function values
+  only once.
+
+  When the reinit function is called, only those gradients, quadrature points
+  etc are transformed to the real cell which belong to the selected face. The
+  number of the selected face is stored such that the #shape_value# function
+  can return the shape function's values on the face which was last selected
+  by a call to the #reinit# function.
+
+  In addition to the complications described above, we need two different
+  Jacobi matrices and determinant in this context: one for the transformation
+  of the unit cell to the real cell (this Jacobi matrix is needed to
+  compute the restriction of the real gradient to the given face) and one
+  for the transformation of the unit face to the real face (needed to
+  compute the weight factors for integration along faces). These two
+  concepts have to be carefully separated.
   */
 template <int dim>
 class FEFaceValues {
@@ -426,7 +426,7 @@ class FEFaceValues {
                                      */
     FEFaceValues (const FiniteElement<dim> &,
                  const Quadrature<dim-1> &,
-                 const UpdateFields);
+                 const UpdateFlags);
 
                                     /**
                                      * Return the value of the #i#th shape
@@ -516,20 +516,28 @@ class FEFaceValues {
                                     /**
                                      * Return the Jacobi determinant times
                                      * the weight of the #i#th quadrature
-                                     * point.
+                                     * point. The Jacobi determinant is that
+                                     * of the transformation of the unit
+                                     * face to the real face, not of the
+                                     * alike cells.
                                      */
     double JxW (const unsigned int i) const;
 
                                     /**
                                      * Return a pointer to the array holding
                                      * the JxW values at the different
-                                     * quadrature points.
+                                     * quadrature points. The Jacobi
+                                     * determinant is that
+                                     * of the transformation of the unit
+                                     * face to the real face, not of the
+                                     * alike cells.
                                      */
     const vector<double> & get_JxW_values () const;
     
                                     /**
                                      * Reinitialize the gradients, Jacobi
-                                     * determinants, etc for the given cell
+                                     * determinants, etc for the face with
+                                     * number #face_no# of #cell#
                                      * and the given finite element.
                                      *
                                      * The constructor needs a boundary object
@@ -545,9 +553,10 @@ class FEFaceValues {
                                      * segments, but higher order elements
                                      * may use other ways.)
                                      */
-    void reinit (const Triangulation<dim>::face_iterator &,
-                const FiniteElement<dim> &,
-                const Boundary<dim> &);
+    void reinit (const DoFHandler<dim>::cell_iterator &cell,
+                const unsigned int                    face_no,
+                const FiniteElement<dim>             &fe,
+                const Boundary<dim>                  &boundary);
 
                                     /**
                                      * Exception
@@ -564,6 +573,14 @@ class FEFaceValues {
                                      * Exception
                                      */
     DeclException0 (ExcCannotInitializeField);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInternalError);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
     
   private:
                                     /**
@@ -573,8 +590,10 @@ class FEFaceValues {
                                      * shape function at the different points,
                                      * columns are for a single point with the
                                      * different shape functions.
+                                     *
+                                     * There is one matrix for each face.
                                      */
-    dFMatrix             shape_values;
+    dFMatrix             shape_values[2*dim];
 
                                     /**
                                      * Store the gradients of the shape
@@ -587,7 +606,11 @@ class FEFaceValues {
                                      * This field is reset each time
                                      * #reinit# is called and contains the
                                      * gradients on the real element, rather
-                                     * than on the reference element.
+                                     * than on the reference element. This
+                                     * function does the transformation from
+                                     * the unit cell to the real cell using
+                                     * the #unit_shape_gradients# for the
+                                     * selected face.
                                      */
     vector<vector<Point<dim> > >  shape_gradients;
 
@@ -598,13 +621,18 @@ class FEFaceValues {
                                      * This field is set up upon construction
                                      * of the object and contains the gradients
                                      * on the reference element.
+                                     *
+                                     * There is one element for each face.
                                      */
-    vector<vector<Point<dim> > >   unit_shape_gradients;
+    vector<vector<Point<dim> > >   unit_shape_gradients[2*dim];
     
                                     /**
                                      * Store an array of the weights of the
                                      * quadrature points. This array is
                                      * set up upon construction.
+                                     *
+                                     * Since these weights are not transformed
+                                     * they are the same for all faces.
                                      */
     vector<double>       weights;
 
@@ -629,13 +657,26 @@ class FEFaceValues {
                                      */
     vector<Point<dim> >  quadrature_points;
 
+                                    /**
+                                     * Array of quadrature points on the
+                                     * unit face. This is a copy of the
+                                     * alike field of the quadrature formula
+                                     * passed upon construction.
+                                     */
+    vector<Point<dim-1> > unit_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.
+                                     *
+                                     * There is one element for each face. The
+                                     * points are computed from those on the
+                                     * unit face, but are stored as coordinates
+                                     * on the unit cell.
                                      */
-    vector<Point<dim-1> > unit_quadrature_points;
+    vector<Point<dim> >  global_unit_quadrature_points[2*dim];
     
                                     /**
                                      * Array of points denoting the off-point
@@ -650,14 +691,35 @@ class FEFaceValues {
                                      * Store the jacobi matrices at the
                                      * different quadrature points. This field
                                      * is set each time #reinit# is called.
+                                     * This is the Jacobi matrix of the
+                                     * transformation of the unit cell to the
+                                     * real cell, not of the unit face to the
+                                     * face. We need this full matrix for the
+                                     * transformation of the gradients to the
+                                     * real cell.
                                      */
     vector<dFMatrix>     jacobi_matrices;
 
+                                    /**
+                                     * List of values denoting the determinant
+                                     * of the transformation from the unit face
+                                     * to the real face. Neede to actually
+                                     * compute the JxW values.
+                                     */
+    vector<double>       face_jacobi_determinants;
+    
                                     /**
                                      * Store which fields are to be updated by
                                      * the reinit function.
                                      */
-    UpdateFields         update_flags;
+    UpdateFlags          update_flags;
+
+                                    /**
+                                     * Store the number of the face selected
+                                     * last time the #reinit# function was
+                                     * called.
+                                     */
+    unsigned int         selected_face;
 };
 
 
@@ -681,7 +743,7 @@ template <int dim>
 inline
 const vector<vector<Point<dim> > > &
 FEValues<dim>::get_shape_grads () const {
-  Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
+  Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
   return shape_gradients;
 };
 
@@ -691,7 +753,7 @@ template <int dim>
 inline
 const vector<Point<dim> > &
 FEValues<dim>::get_quadrature_points () const {
-  Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
+  Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   return quadrature_points;
 };
 
@@ -701,7 +763,7 @@ template <int dim>
 inline
 const vector<Point<dim> > &
 FEValues<dim>::get_ansatz_points () const {
-  Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
+  Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   return ansatz_points;
 };
 
@@ -711,7 +773,7 @@ template <int dim>
 inline
 const vector<double> &
 FEValues<dim>::get_JxW_values () const {
-  Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
+  Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   return JxW_values;
 };
 
@@ -725,7 +787,7 @@ FEValues<dim>::get_JxW_values () const {
 template <int dim>
 inline
 const dFMatrix & FEFaceValues<dim>::get_shape_values () const {
-  return shape_values;
+  return shape_values[selected_face];
 };
 
 
@@ -735,7 +797,7 @@ template <int dim>
 inline
 const vector<vector<Point<dim> > > &
 FEFaceValues<dim>::get_shape_grads () const {
-  Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
+  Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
   return shape_gradients;
 };
 
@@ -745,7 +807,7 @@ template <int dim>
 inline
 const vector<Point<dim> > &
 FEFaceValues<dim>::get_quadrature_points () const {
-  Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
+  Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   return quadrature_points;
 };
 
@@ -755,7 +817,7 @@ template <int dim>
 inline
 const vector<Point<dim> > &
 FEFaceValues<dim>::get_ansatz_points () const {
-  Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
+  Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   return ansatz_points;
 };
 
@@ -765,7 +827,7 @@ template <int dim>
 inline
 const vector<double> &
 FEFaceValues<dim>::get_JxW_values () const {
-  Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
+  Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   return JxW_values;
 };
 
index f4a962c8472720f4a3777be67ef3f156f0ebf40d..e08fda1ac9576a5e5fdcc98b8d2bc020866acfea 100644 (file)
@@ -3,26 +3,76 @@
 #include <fe/fe.h>
 #include <fe/quadrature.h>
 #include <grid/tria_iterator.h>
-#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
 #include <grid/tria_boundary.h>
 
 
 
 
 
+/*------------------------------- FiniteElementData ----------------------*/
+
+
+bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const {
+  return ((dofs_per_vertex == f.dofs_per_vertex) &&
+         (dofs_per_line == f.dofs_per_line) &&
+         (total_dofs == f.total_dofs));
+};
+
+
+
+bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const {
+  return ((dofs_per_vertex == f.dofs_per_vertex) &&
+         (dofs_per_line == f.dofs_per_line) &&
+         (dofs_per_quad == f.dofs_per_quad) &&
+         (total_dofs == f.total_dofs));
+};
+
+
+
+
+
 /*------------------------------- FiniteElementBase ----------------------*/
 
-template <int dim>
-FiniteElementBase<dim>::FiniteElementBase (const FiniteElementBase<dim> &f) :
-               total_dofs(f.total_dofs),
-               interface_constraints (f.interface_constraints)
+
+
+
+FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex,
+                                        const unsigned int dofs_per_line,
+                                        const unsigned int dofs_per_quad) :
+               FiniteElementData<1> (dofs_per_vertex,
+                                     dofs_per_line)
+{
+  Assert (dofs_per_quad==0, ExcInternalError());
+
+  const unsigned int dim=1;
+  for (unsigned int i=0; i<(1<<dim); ++i) 
+    {
+      restriction[i].reinit (total_dofs, total_dofs);
+      prolongation[i].reinit (total_dofs, total_dofs);
+    };
+  interface_constraints.reinit (1,1);
+  interface_constraints(0,0)=1.;
+};
+
+
+
+FiniteElementBase<2>::FiniteElementBase (const unsigned int dofs_per_vertex,
+                                        const unsigned int dofs_per_line,
+                                        const unsigned int dofs_per_quad) :
+               FiniteElementData<2> (dofs_per_vertex,
+                                     dofs_per_line,
+                                     dofs_per_quad)
 {
+  const unsigned int dim=2;
   for (unsigned int i=0; i<(1<<dim); ++i) 
     {
-      restriction[i] = f.restriction[i];
-      prolongation[i] = f.prolongation[i];
+      restriction[i].reinit (total_dofs, total_dofs);
+      prolongation[i].reinit (total_dofs, total_dofs);
     };
-}
+  interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line,
+                               2*dofs_per_vertex+dofs_per_line);
+};
 
 
 
@@ -47,6 +97,11 @@ FiniteElementBase<dim>::prolongate (const unsigned int child) const {
 template <int dim>
 const dFMatrix &
 FiniteElementBase<dim>::constraints () const {
+  if (dim==1)
+    Assert ((interface_constraints.m()==1) && (interface_constraints.n()==1),
+           ExcWrongInterfaceMatrixSize(interface_constraints.m(),
+                                       interface_constraints.n()));
+  
   return interface_constraints;
 };
 
@@ -54,67 +109,19 @@ FiniteElementBase<dim>::constraints () const {
 
 template <int dim>
 bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const {
-  return ((total_dofs == f.total_dofs) &&
+  return ((static_cast<FiniteElementData<dim> >(*this) ==
+          static_cast<FiniteElementData<dim> >(f)) &&
          (interface_constraints == f.interface_constraints));
 };
 
 
 
-template <int dim>
-double FiniteElementBase<dim>::shape_value (const unsigned int,
-                                           const Point<dim> &) const {
-  Assert (false, ExcPureFunctionCalled());
-  return 0;
-};
-
-
-
-template <int dim>
-Point<dim> FiniteElementBase<dim>::shape_grad (const unsigned int,
-                                              const Point<dim> &) const {
-  Assert (false, ExcPureFunctionCalled());
-  return Point<dim>();
-};
-
-
-
-
-template <int dim>
-void FiniteElementBase<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator &,
-                                            const vector<Point<dim> > &,
-                                            vector<dFMatrix> &,
-                                            const bool,
-                                            vector<Point<dim> > &,
-                                            const bool,
-                                            vector<Point<dim> > &,
-                                            const bool,
-                                            const Boundary<dim> &) const {
-  Assert (false, ExcPureFunctionCalled());
-};
-
-
-
-template <int dim>
-void FiniteElementBase<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &,
-                                                const Boundary<dim> &,
-                                                vector<Point<dim> > &) const {
-  Assert (false, ExcPureFunctionCalled());
-};
-  
 
 
 /*------------------------------- FiniteElement ----------------------*/
 
 
-bool FiniteElement<1>::operator == (const FiniteElement<1> &f) const {
-  return ((dofs_per_vertex == f.dofs_per_vertex) &&
-         (dofs_per_line == f.dofs_per_line) &&
-         (FiniteElementBase<1>::operator == (f)));
-};
-
-
-
-void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
+void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
                                       const vector<Point<1> > &unit_points,
                                       vector<dFMatrix>  &jacobians,
                                       const bool         compute_jacobians,
@@ -162,54 +169,97 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce
 };
 
 
-
-void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &,
-                                          const Boundary<1> &,
-                                          vector<Point<1> > &) const {
-                                  // is this function useful in 1D?
+template <int dim>
+void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
+                                        const vector<Point<dim> > &,
+                                        vector<dFMatrix> &,
+                                        const bool,
+                                        vector<Point<dim> > &,
+                                        const bool,
+                                        vector<Point<dim> > &,
+                                        const bool,
+                                        const Boundary<dim> &) const {
   Assert (false, ExcPureFunctionCalled());
 };
 
 
 
-bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const {
-  return ((dofs_per_vertex == f.dofs_per_vertex) &&
-         (dofs_per_line == f.dofs_per_line) &&
-         (dofs_per_quad == f.dofs_per_quad) &&
-         (FiniteElementBase<2>::operator == (f)));
-};
-
+void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &,
+                                           const unsigned int       ,
+                                           const vector<Point<0> > &,
+                                           const vector<Point<1> > &,
+                                           vector<dFMatrix>        &,
+                                           const bool               ,
+                                           vector<Point<1> >       &,
+                                           const bool               ,
+                                           vector<Point<1> >       &,
+                                           const bool               ,
+                                           vector<double>          &,
+                                           const bool              ,
+                                           const Boundary<1>       &) const {
+  Assert (false, ExcNotImplemented());
+}
 
 
-void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
-                                      const vector<Point<2> > &,
-                                      vector<dFMatrix> &,
-                                      const bool,
-                                      vector<Point<2> > &,
-                                      const bool,
-                                      vector<Point<2> > &,
-                                      const bool,
-                                      const Boundary<2> &) const {
-  Assert (false, ExcPureFunctionCalled());
-};
 
+template <int dim>
+void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterator &cell,
+                                             const unsigned int           face_no,
+                                             const vector<Point<dim-1> > &unit_points,
+                                             const vector<Point<dim> > &global_unit_points,
+                                             vector<dFMatrix>    &jacobians,
+                                             const bool           compute_jacobians,
+                                             vector<Point<dim> > &ansatz_points,
+                                             const bool           compute_ansatz_points,
+                                             vector<Point<dim> > &q_points,
+                                             const bool           compute_q_points,
+                                             vector<double>      &face_jacobi_determinants,
+                                             const bool           compute_face_jacobians,
+                                             const Boundary<dim> &boundary) const {
+  Assert (jacobians.size() == unit_points.size(),
+         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+  Assert (q_points.size() == unit_points.size(),
+         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+  Assert (global_unit_points.size() == unit_points.size(),
+         ExcWrongFieldDimension(global_unit_points.size(), unit_points.size()));
+  Assert (ansatz_points.size() == dofs_per_face,
+         ExcWrongFieldDimension(ansatz_points.size(), dofs_per_face));
+  
+  static vector<Point<dim> > dummy(total_dofs);
+  fill_fe_values (cell, global_unit_points,
+                 jacobians, compute_jacobians,
+                 dummy, false,
+                 q_points, compute_q_points,
+                 boundary);
+  
+  cout << "Global unit points:\n";
+  for (unsigned int p=0; p<unit_points.size(); ++p)
+    cout << "    " << global_unit_points[p] << endl;
 
+  if (compute_ansatz_points)
+    get_face_ansatz_points (cell->face(face_no), boundary, ansatz_points);
 
-void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &,
-                                          const Boundary<2> &,
-                                          vector<Point<2> > &) const {
-                                  // is this function useful in 1D?
-  Assert (false, ExcPureFunctionCalled());
+  if (compute_face_jacobians)
+    get_face_jacobians (cell->face(face_no), boundary,
+                       unit_points, face_jacobi_determinants);
+  
+  cout << "Global ansatz points:\n";
+  for (unsigned int p=0; p<unit_points.size(); ++p)
+    cout << "    " << ansatz_points[p] << endl;
 };
 
 
 
-
-
 /*------------------------------- Explicit Instantiations -------------*/
 
+template class FiniteElementData<1>;
+template class FiniteElementData<2>;
+
 template class FiniteElementBase<1>;
 template class FiniteElementBase<2>;
 
 template class FiniteElement<1>;
 template class FiniteElement<2>;
+
+
+
index 5eafa7b819963012ed39ba43f64b35cbcf9a6bdd..45313af52c2abb58de76c80561510f78390c188a 100644 (file)
@@ -2,17 +2,14 @@
 
 #include <fe/fe_lib.h>
 #include <grid/tria_iterator.h>
-#include <grid/tria_accessor.h>
-
+#include <grid/dof_accessor.h>
+#include <algorithm>
 
 
 
 FELinear<1>::FELinear () :
                FiniteElement<1> (1, 0)
 {
-  restriction[0].reinit (2,2);
-  restriction[1].reinit (2,2);
-
                                   // for restriction and prolongation matrices:
                                   // note that we do not add up all the
                                   // contributions since then we would get
@@ -33,10 +30,6 @@ FELinear<1>::FELinear () :
   restriction[1](1,0) = 1./2.;
   restriction[1](1,1) = 1.0;
 
-
-  prolongation[0].reinit (2,2);
-  prolongation[1].reinit (2,2);
-  
   prolongation[0](0,0) = 1.0;
   prolongation[0](1,0) = 1./2.;
   prolongation[0](1,1) = 1./2.;
@@ -78,8 +71,8 @@ FELinear<1>::shape_grad(const unsigned int i,
 
 
 
-void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
-                                 const vector<Point<1> >               &unit_points,
+void FELinear<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+                                 const vector<Point<1> >            &unit_points,
                                  vector<dFMatrix>  &jacobians,
                                  const bool         compute_jacobians,
                                  vector<Point<1> > &ansatz_points,
@@ -96,32 +89,29 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
 
 
 
-void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &,
-                                     const Boundary<1> &,
-                                     vector<Point<1> > &) const {
-  Assert (false, ExcPureFunctionCalled());
+void FELinear<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+                                         const Boundary<1>  &,
+                                         vector<Point<1> >  &) const {
+  Assert (false, ExcInternalError());
 };
 
 
 
+void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+                                     const Boundary<1>         &,
+                                     const vector<Point<0> > &,
+                                     vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
 
 FELinear<2>::FELinear () :
                FiniteElement<2> (1, 0, 0)
 {
-  interface_constraints.reinit(1,2);
   interface_constraints(0,0) = 1./2.;
   interface_constraints(0,1) = 1./2.;
 
-  restriction[0].reinit(4,4);
-  restriction[1].reinit(4,4);
-  restriction[2].reinit(4,4);
-  restriction[3].reinit(4,4);
-
-  prolongation[0].reinit(4,4);
-  prolongation[1].reinit(4,4);
-  prolongation[2].reinit(4,4);
-  prolongation[3].reinit(4,4);
-
   restriction[0](0,0) = 1.0;
   restriction[0](0,1) = 1./2.;
   restriction[0](1,1) = 1./2.;
@@ -239,17 +229,16 @@ FELinear<2>::shape_grad (const unsigned int i,
 
 
 
-// this function may be generalised to three or more dimensions with gcc2.8
-// you will have to change the number of vertices
-void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
-                                 const vector<Point<2> >               &unit_points,
-                                 vector<dFMatrix>  &jacobians,
-                                 const bool         compute_jacobians,
-                                 vector<Point<2> > &ansatz_points,
-                                 const bool         compute_ansatz_points,
-                                 vector<Point<2> > &q_points,
-                                 const bool         compute_q_points,
-                                 const Boundary<2> &) const {
+template <int dim>
+void FELinear<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                   const vector<Point<dim> >            &unit_points,
+                                   vector<dFMatrix>    &jacobians,
+                                   const bool           compute_jacobians,
+                                   vector<Point<dim> > &ansatz_points,
+                                   const bool           compute_ansatz_points,
+                                   vector<Point<dim> > &q_points,
+                                   const bool           compute_q_points,
+                                   const Boundary<dim> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -257,7 +246,6 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
   Assert (ansatz_points.size() == total_dofs,
          ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
   
-  const unsigned int dim=2;
   const unsigned int n_vertices=(1<<dim);
   unsigned int n_points=unit_points.size();
 
@@ -330,26 +318,44 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                   // compute ansatz points, which are
                                   // the corners for linear elements
   if (compute_ansatz_points) 
-    for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int vertex=0; vertex<n_vertices; ++vertex)
       ansatz_points[vertex] = vertices[vertex];
 };
 
 
 
+
 template <int dim>
-void FELinear<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &face,
-                                       const Boundary<dim>  &,
-                                       vector<Point<dim> >  &ansatz_points) const {
+void FELinear<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face_iterator &face,
+                                           const Boundary<dim>  &,
+                                           vector<Point<dim> >  &ansatz_points) const {
   Assert (ansatz_points.size() == (1<<(dim-1)),
-         typename FiniteElementBase<dim>::ExcWrongFieldDimension (ansatz_points.size(),
-                                                                  1<<(dim-1)));
-  
+         ExcWrongFieldDimension (ansatz_points.size(), 1<<(dim-1)));
+
   for (unsigned int vertex=0; vertex<(1<<(dim-1)); ++vertex)
     ansatz_points[vertex] = face->vertex(vertex);
 };
 
 
 
+void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
+                                     const Boundary<2>         &,
+                                     const vector<Point<1> > &unit_points,
+                                     vector<double>      &face_jacobians) const {
+  Assert (unit_points.size() == face_jacobians.size(),
+         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
+
+                                  // a linear mapping for a single line
+                                  // produces particularly simple
+                                  // expressions for the jacobi
+                                  // determinant :-)
+  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
+  fill_n (face_jacobians.begin(),
+         unit_points.size(),
+         h);  
+};
+
+
 
 
 
@@ -359,8 +365,8 @@ FEQuadratic<1>::FEQuadratic () :
 
 
 
-void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
-                                    const vector<Point<1> >               &unit_points,
+void FEQuadratic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+                                    const vector<Point<1> >            &unit_points,
                                     vector<dFMatrix>  &jacobians,
                                     const bool         compute_jacobians,
                                     vector<Point<1> > &ansatz_points,
@@ -377,10 +383,26 @@ void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell
 
 
 
+void FEQuadratic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+                                            const Boundary<1>  &,
+                                            vector<Point<1> >  &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+                                        const Boundary<1>         &,
+                                        const vector<Point<0> > &,
+                                        vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 FEQuadratic<2>::FEQuadratic () :
                FiniteElement<2> (1, 1, 1)
 {
-  interface_constraints.reinit(3,3);
   interface_constraints(0,2) = 1.0;
   interface_constraints(1,0) = 3./8.;
   interface_constraints(1,1) = -1./8.;
@@ -388,6 +410,10 @@ FEQuadratic<2>::FEQuadratic () :
   interface_constraints(2,0) = -1./8.;
   interface_constraints(2,1) = 3./8.;
   interface_constraints(2,2) = 3./4.;
+
+                                  // still implement restriction
+                                  // and prolongation
+  Assert (false, ExcNotImplemented());
 };
 
 
@@ -398,7 +424,7 @@ FEQuadratic<dim>::shape_value (const unsigned int i,
                               const Point<dim> &) const
 {
   Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
   return 0.;
 };
 
@@ -410,21 +436,22 @@ FEQuadratic<dim>::shape_grad (const unsigned int i,
                              const Point<dim> &) const
 {
   Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
   return Point<dim> ();
 };
 
 
 
-void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
-                                    const vector<Point<2> >               &unit_points,
-                                    vector<dFMatrix>  &jacobians,
-                                    const bool,
-                                    vector<Point<2> > &ansatz_points,
-                                    const bool,
-                                    vector<Point<2> > &q_points,
-                                    const bool,
-                                    const Boundary<2> &) const {
+template <int dim>
+void FEQuadratic<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
+                                      const vector<Point<dim> >            &unit_points,
+                                      vector<dFMatrix>  &jacobians,
+                                      const bool,
+                                      vector<Point<dim> > &ansatz_points,
+                                      const bool,
+                                      vector<Point<dim> > &q_points,
+                                      const bool,
+                                      const Boundary<dim> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -432,7 +459,26 @@ void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
   Assert (ansatz_points.size() == total_dofs,
          ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
 
-  Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face_iterator &,
+                                              const Boundary<dim>  &,
+                                              vector<Point<dim> >  &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
+                                     const Boundary<dim>         &,
+                                     const vector<Point<dim-1> > &,
+                                     vector<double>      &) const {
+  Assert (false, ExcNotImplemented());
 };
 
 
@@ -445,8 +491,8 @@ FECubic<1>::FECubic () :
 
 
 
-void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
-                                const vector<Point<1> >               &unit_points,
+void FECubic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+                                const vector<Point<1> >            &unit_points,
                                 vector<dFMatrix>  &jacobians,
                                 const bool         compute_jacobians,
                                 vector<Point<1> > &ansatz_points,
@@ -463,6 +509,23 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
 
 
 
+void FECubic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+                                        const Boundary<1>  &,
+                                        vector<Point<1> >  &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+                                    const Boundary<1>         &,
+                                    const vector<Point<0> > &,
+                                    vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 
 FECubic<2>::FECubic () :
                FiniteElement<2> (1, 2, 4) {};
@@ -476,7 +539,7 @@ FECubic<dim>::shape_value (const unsigned int i,
                           const Point<dim> &) const
 {
   Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
   return 0.;
 };
 
@@ -488,21 +551,22 @@ FECubic<dim>::shape_grad (const unsigned int i,
                          const Point<dim> &) const
 {
   Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
   return Point<dim> ();
 };
 
 
 
-void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
-                                const vector<Point<2> >               &unit_points,
-                                vector<dFMatrix>  &jacobians,
-                                const bool,
-                                vector<Point<2> > &ansatz_points,
-                                const bool,
-                                vector<Point<2> > &q_points,
-                                const bool,
-                                const Boundary<2> &) const {
+template <int dim>
+void FECubic<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
+                                  const vector<Point<dim> >            &unit_points,
+                                  vector<dFMatrix>  &jacobians,
+                                  const bool,
+                                  vector<Point<dim> > &ansatz_points,
+                                  const bool,
+                                  vector<Point<dim> > &q_points,
+                                  const bool,
+                                  const Boundary<dim> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -510,15 +574,32 @@ void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
   Assert (ansatz_points.size() == total_dofs,
          ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
 
-  Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
+  Assert (false, ExcNotImplemented());
 };
 
 
 
+template <int dim>
+void FECubic<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face_iterator &,
+                                          const Boundary<dim>  &,
+                                          vector<Point<dim> >  &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FECubic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
+                                      const Boundary<dim>         &,
+                                      const vector<Point<dim-1> > &,
+                                      vector<double>      &) const {
+  Assert (false, ExcNotImplemented());
+};
+
 
 
 
-// explicite instantiations
+// explicit instantiations
 
 template class FELinear<1>;
 template class FELinear<2>;
index dae361ba7a463befc62574beb1f37554fb124922..a332c178a432f47cd2e8d3d209b2dce9eba2d6c5 100644 (file)
@@ -15,7 +15,7 @@
 template <int dim>
 FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                         const Quadrature<dim>    &quadrature,
-                        const UpdateFields        update_flags) :
+                        const UpdateFlags         update_flags) :
                n_quadrature_points(quadrature.n_quadrature_points),
                total_dofs(fe.total_dofs),
                shape_values(fe.total_dofs, quadrature.n_quadrature_points),
@@ -66,7 +66,7 @@ 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());
+  Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
 
   return shape_gradients[i][j];
 };
@@ -76,7 +76,7 @@ FEValues<dim>::shape_grad (const unsigned int 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));
-  Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
+  Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   
   return quadrature_points[i];
 };
@@ -86,7 +86,7 @@ const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const
 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());
+  Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   
   return ansatz_points[i];
 };
@@ -96,7 +96,7 @@ const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
 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());
+  Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   
   return JxW_values[i];
 };
@@ -104,28 +104,29 @@ double FEValues<dim>::JxW (const unsigned int i) const {
 
 
 template <int dim>
-void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
-                           const FiniteElement<dim>                         &fe,
-                           const Boundary<dim>                              &boundary) {
+void FEValues<dim>::reinit (const typename DoFHandler<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))
+  if ((update_flags & update_jacobians) ||
+      (update_flags & update_q_points)  ||
+      (update_flags & update_ansatz_points))
     fe.fill_fe_values (cell,
                       unit_quadrature_points,
                       jacobi_matrices,
-                      update_flags | update_jacobians,
+                      update_flags & update_jacobians,
                       ansatz_points,
-                      update_flags | update_ansatz_points,
+                      update_flags & update_ansatz_points,
                       quadrature_points,
-                      update_flags | update_q_points,
+                      update_flags & update_q_points,
                       boundary);
 
                                   // compute gradients on real element if
                                   // requested
-  if (update_flags | update_gradients) 
+  if (update_flags & update_gradients) 
     {
-      Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+      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)
@@ -149,9 +150,9 @@ void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &ce
                                   // refer to the general doc for
                                   // why we take the inverse of the
                                   // determinant
-  if (update_flags | update_JxW_values) 
+  if (update_flags & update_JxW_values) 
     {
-      Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+      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();
     };
@@ -167,36 +168,84 @@ void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &ce
 template <int dim>
 FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                                 const Quadrature<dim-1>  &quadrature,
-                                const UpdateFields        update_flags) :
+                                const UpdateFlags         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-1>()),
                ansatz_points (fe.total_dofs, Point<dim>()),
-               jacobi_matrices (quadrature.n_quadrature_points,
-                                dFMatrix(dim,dim)),
-               update_flags (update_flags)
+               jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)),
+               face_jacobi_determinants (quadrature.n_quadrature_points,0),
+               update_flags (update_flags),
+               selected_face(0)
 {
-  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 face=0; face<2*dim; ++face)
+    {
+      shape_values[face].reinit(fe.total_dofs, quadrature.n_quadrature_points);
+      unit_shape_gradients[face].resize (fe.total_dofs,
+                                        vector<Point<dim> >(quadrature.n_quadrature_points));
+      global_unit_quadrature_points[face].resize (quadrature.n_quadrature_points,
+                                                 Point<dim>());
+    };
+
+                                  // set up an array of the unit points
+                                  // on the given face, but in coordinates
+                                  // of the space with #dim# dimensions.
+                                  // the points are still on the unit
+                                  // cell.
+  for (unsigned int face=0; face<2*dim; ++face)
+    for (unsigned int p=0; p<n_quadrature_points; ++p)
+      switch (dim) 
+       {
+         case 2:
+         {
+           
+           switch (face)
+             {
+               case 0:
+                     global_unit_quadrature_points[face][p]
+                       = Point<dim>(unit_quadrature_points[p](0),0);
+                     break;       
+               case 1:
+                     global_unit_quadrature_points[face][p]
+                       = Point<dim>(1,unit_quadrature_points[p](0));
+                     break;       
+               case 2:
+                     global_unit_quadrature_points[face][p]
+                       = Point<dim>(unit_quadrature_points[p](0),1);
+                     break;       
+               case 3:
+                     global_unit_quadrature_points[face][p]
+                       = Point<dim>(0,unit_quadrature_points[p](0));
+                     break;
+               default:
+                     Assert (false, ExcInternalError());
+             };
+           
+           break;
+         };
+         default:
+               Assert (false, ExcNotImplemented());
+       };
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
       weights[i] = quadrature.weight(i);
       unit_quadrature_points[i] = quadrature.quad_point(i);
     };
+
+  for (unsigned int face=0; face<2*dim; ++face)
+    for (unsigned int i=0; i<fe.total_dofs; ++i)
+      for (unsigned int j=0; j<n_quadrature_points; ++j) 
+       {
+         shape_values[face](i,j) = fe.shape_value(i, global_unit_quadrature_points[face][j]);
+         unit_shape_gradients[face][i][j]
+           = fe.shape_grad(i, global_unit_quadrature_points[face][j]);
+       };
 };
 
 
@@ -204,10 +253,12 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
 template <int dim>
 double FEFaceValues<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()));
+  Assert (i<shape_values[selected_face].m(),
+         ExcInvalidIndex (i, shape_values[selected_face].m()));
+  Assert (j<shape_values[selected_face].n(),
+         ExcInvalidIndex (j, shape_values[selected_face].n()));
 
-  return shape_values(i,j);
+  return shape_values[selected_face](i,j);
 };
 
 
@@ -216,9 +267,11 @@ template <int dim>
 const Point<dim> &
 FEFaceValues<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());
+  Assert (i<shape_values[selected_face].m(),
+         ExcInvalidIndex (i, shape_values[selected_face].m()));
+  Assert (j<shape_values[selected_face].n(),
+         ExcInvalidIndex (j, shape_values[selected_face].n()));
+  Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
 
   return shape_gradients[i][j];
 };
@@ -228,7 +281,7 @@ FEFaceValues<dim>::shape_grad (const unsigned int i,
 template <int dim>
 const Point<dim> & FEFaceValues<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());
+  Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   
   return quadrature_points[i];
 };
@@ -238,7 +291,7 @@ const Point<dim> & FEFaceValues<dim>::quadrature_point (const unsigned int i) co
 template <int dim>
 const Point<dim> & FEFaceValues<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());
+  Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   
   return ansatz_points[i];
 };
@@ -248,7 +301,7 @@ const Point<dim> & FEFaceValues<dim>::ansatz_point (const unsigned int i) const
 template <int dim>
 double FEFaceValues<dim>::JxW (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
-  Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
+  Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   
   return JxW_values[i];
 };
@@ -256,29 +309,36 @@ double FEFaceValues<dim>::JxW (const unsigned int i) const {
 
 
 template <int dim>
-void FEFaceValues<dim>::reinit (const typename Triangulation<dim>::face_iterator &face,
-                               const FiniteElement<dim>                         &fe,
-                               const Boundary<dim>                              &boundary) {
+void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
+                               const unsigned int                             face_no,
+                               const FiniteElement<dim>                      &fe,
+                               const Boundary<dim>                           &boundary) {
+  selected_face = face_no;
                                   // fill jacobi matrices and real
                                   // quadrature points
-  if ((update_flags | update_jacobians) ||
-      (update_flags | update_q_points))
-//**    fe.fill_fe_values (face,
-//**                  unit_quadrature_points,
-//**                  jacobi_matrices,
-//**                  update_flags | update_jacobians,
-//**                  ansatz_points,
-//**                  update_flags | update_ansatz_points,
-//**                  quadrature_points,
-//**                  update_flags | update_q_points,
-//**                  boundary)
-      ;
+  if ((update_flags & update_jacobians) ||
+      (update_flags & update_q_points)  ||
+      (update_flags & update_ansatz_points) ||
+      (update_flags & update_JxW_values))
+    fe.fill_fe_face_values (cell,
+                           face_no,
+                           unit_quadrature_points,
+                           global_unit_quadrature_points[face_no],
+                           jacobi_matrices,
+                           update_flags & update_jacobians,
+                           ansatz_points,
+                           update_flags & update_ansatz_points,
+                           quadrature_points,
+                           update_flags & update_q_points,
+                           face_jacobi_determinants,
+                           update_flags & update_JxW_values,
+                           boundary);
 
                                   // compute gradients on real element if
                                   // requested
-  if (update_flags | update_gradients) 
+  if (update_flags & update_gradients) 
     {
-      Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+      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)
@@ -292,7 +352,7 @@ void FEFaceValues<dim>::reinit (const typename Triangulation<dim>::face_iterator
              for (unsigned int b=0; b<dim; ++b)
                shape_gradients[i][j](s)
                  +=
-                 unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
+                 unit_shape_gradients[face_no][i][j](b) * jacobi_matrices[j](b,s);
            };
     };
   
@@ -302,11 +362,11 @@ void FEFaceValues<dim>::reinit (const typename Triangulation<dim>::face_iterator
                                   // refer to the general doc for
                                   // why we take the inverse of the
                                   // determinant
-  if (update_flags | update_JxW_values) 
+  if (update_flags & update_JxW_values) 
     {
-      Assert (update_flags | update_jacobians, ExcCannotInitializeField());
+      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();
+       JxW_values[i] = weights[i] * face_jacobi_determinants[i];
     };
 };
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.