]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add more support for curved boundaries and the integration on faces.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 24 Apr 1998 14:13:33 +0000 (14:13 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 24 Apr 1998 14:13:33 +0000 (14:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@191 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Todo
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/include/numerics/assembler.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc

index 859ab7724eb871a86f80923c173bd9dcf4eea6ac..a235c34cee9b2f4d4bfbbacbb408ae19cfc953cc 100644 (file)
@@ -20,12 +20,17 @@ Write monitors to control whether enough memory was allocated for
 
 Write Triangulation::mesh_function.
 
+
 Unify CellAccessor<1> and <2> by renaming
     LineAccessor<dim> -> SubstructAccessor<dim,1>
     QuadAccessor<dim> -> SubstructAccessor<dim,2>
   and deriving CellAccessor<dim> from SubstructAccessor<dim,dim>
   Do the same with DoFLineAccessor and DoFQuadAccessor
 
+Unify FEValues<dim> and FEFaceValues<dim> to FEValues<dim,subdim>
+  and use partial specialization.
+
+
 Add support for fast assemblage of mass matrices.
 
 Make AssemblerData a local class to Assembler again if gcc2.8 supports it.
index 482e694caab77d3eaeebaaa62294c9c117200c71..8a81163cbbb1ce9ca387b915070cbc2a1ac737dc 100644 (file)
@@ -12,6 +12,7 @@
 
 
 // forward declarations
+template <int dim> class Boundary;
 template <int dim> class FiniteElement;
 template <int dim> class Quadrature;
 
@@ -127,6 +128,13 @@ class FEValues {
                                      * specified finite element using the
                                      * quadrature points of the given
                                      * quadrature rule.
+                                     *
+                                     * This function actually only fills
+                                     * the fields related to the unit face,
+                                     * the fields related to a real face (like
+                                     * gradients, true quadrature points, etc.)
+                                     * need to be initialized using the
+                                     * #reinit# function.
                                      */
     FEValues (const FiniteElement<dim> &,
              const Quadrature<dim> &,
@@ -228,9 +236,23 @@ class FEValues {
                                      * Reinitialize the gradients, Jacobi
                                      * determinants, etc for the given cell
                                      * and the given finite element.
+                                     *
+                                     * This function needs a boundary object
+                                     * passed, since this class needs to know
+                                     * how to handle faces which are located
+                                     * on the boundary of the domain. In that
+                                     * case, faces may be curved and the
+                                     * calculation of quadrature points,
+                                     * gradients and the like may need
+                                     * additional effort, depending on the
+                                     * mapping from the unit to the real cell
+                                     * (linear mappings use straight boundary
+                                     * segments, but higher order elements
+                                     * may use other ways.)
                                      */
     void reinit (const Triangulation<dim>::cell_iterator &,
-                const FiniteElement<dim> &);
+                const FiniteElement<dim> &,
+                const Boundary<dim> &);
 
                                     /**
                                      * Exception
@@ -346,6 +368,306 @@ class FEValues {
 
 
 
+/**
+  Represent a finite element evaluated with a specific quadrature rule.
+  This class is an optimization which avoids evaluating the shape functions
+  at the quadrature points each time a quadrature takes place. Rather, the
+  values and gradients (and possibly higher order derivatives in future
+  versions of this library) are evaluated once and for all on the unit
+  face before doing the quadrature itself. Only the Jacobian matrix of
+  the transformation from the unit face to the real face and the integration
+  points in real space are calculated each time we move on to a new face.
+
+  The unit face is defined to be the tensor product of the interval $[0,1]$
+  in the present number of dimensions minus one. In part of the literature,
+  the convention is used that the unit cell be the tensor product of the
+  interval $[-1,1]$, which is to distinguished properly.
+
+  This class is very similar to the #FEValues# class; see there for more
+  documentation.
+  */
+template <int dim>
+class FEFaceValues {
+  public:
+
+    
+    
+                                    /**
+                                     * Number of quadrature points on
+                                     * the face.
+                                     */
+    const unsigned int n_quadrature_points;
+
+                                    /**
+                                     * Total number of shape functions
+                                     * on the cell adjacent to this face.
+                                     * This number is not the same as the
+                                     * number of shape functions of which
+                                     * the center is located on the face.
+                                     */
+    const unsigned int total_dofs;
+    
+                                    /**
+                                     * Constructor. Fill all arrays with the
+                                     * values of the shape functions of the
+                                     * specified finite element using the
+                                     * quadrature points of the given
+                                     * quadrature rule for the face, which
+                                     * has a dimension one less than the
+                                     * cell.
+                                     *
+                                     * This function actually only fills
+                                     * the fields related to the unit face,
+                                     * the fields related to a real face (like
+                                     * gradients, true quadrature points, etc.)
+                                     * need to be initialized using the
+                                     * #reinit# function.
+                                     */
+    FEFaceValues (const FiniteElement<dim> &,
+                 const Quadrature<dim-1> &,
+                 const UpdateFields);
+
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at the #j# quadrature point.
+                                     */
+    double shape_value (const unsigned int i,
+                       const unsigned int j) const;
+
+                                    /**
+                                     * Return a pointer to the matrix holding
+                                     * all values of shape functions at all
+                                     * integration points, on the present cell.
+                                     * For the format of this matrix, see the
+                                     * documentation for the matrix itself.
+                                     */
+    const dFMatrix & get_shape_values () const;
+    
+                                    /**
+                                     * Return the gradient of the #i#th shape
+                                     * function at the #j# quadrature point.
+                                     * If you want to get the derivative in
+                                     * one of the coordinate directions, use
+                                     * the appropriate function of the #Point#
+                                     * class to extract one component. Since
+                                     * only a reference to the gradient's value
+                                     * is returned, there should be no major
+                                     * performance drawback.
+                                     * The function returns the gradient on the
+                                     * real element, not the reference element.
+                                     */
+    const Point<dim> & shape_grad (const unsigned int i,
+                                  const unsigned int j) const;
+
+                                    /** 
+                                     * Return a pointer to the matrix holding
+                                     * all gradients of shape functions at all
+                                     * integration points, on the present cell.
+                                     * For the format of this matrix, see the
+                                     * documentation for the matrix itself.
+                                     */
+    const vector<vector<Point<dim> > > & get_shape_grads () const;
+    
+                                    /**
+                                     * Return the position of the #i#th
+                                     * quadrature point in real space.
+                                     *
+                                     * For curved boundary cells, using
+                                     * biquadratic or higher mappings
+                                     * of the unit cell to the real cell,
+                                     * these points may not be on the
+                                     * plane submannifold on which the
+                                     * vertices of the face lie.
+                                     */
+    const Point<dim> & quadrature_point (const unsigned int i) const;
+
+                                    /**
+                                     * Return a pointer to the vector of
+                                     * quadrature points.
+                                     */
+    const vector<Point<dim> > & get_quadrature_points () const;
+
+                                    /**
+                                     * Return the point in real space where
+                                     * the #i#th ansatz function is located
+                                     * (location is in the sense of where it
+                                     * assumes its nominal properties, e.g. at
+                                     * the vertex of a cell, at the center of
+                                     * a line, etc).
+                                     *
+                                     * This function is needed for the
+                                     * interpolation problem: if we want to
+                                     * transfer a continuous function to a
+                                     * finite element function by interpolation
+                                     * we have to take the continuous
+                                     * function's value at the ansatz function
+                                     * locations.
+                                     */
+    const Point<dim> & ansatz_point (const unsigned int i) const;
+
+                                    /**
+                                     * Return a pointer to the vector of points
+                                     * denoting the location of the ansatz
+                                     * functions.
+                                     */
+    const vector<Point<dim> > & get_ansatz_points () const;
+    
+                                    /**
+                                     * Return the Jacobi determinant times
+                                     * the weight of the #i#th quadrature
+                                     * point.
+                                     */
+    double JxW (const unsigned int i) const;
+
+                                    /**
+                                     * Return a pointer to the array holding
+                                     * the JxW values at the different
+                                     * quadrature points.
+                                     */
+    const vector<double> & get_JxW_values () const;
+    
+                                    /**
+                                     * Reinitialize the gradients, Jacobi
+                                     * determinants, etc for the given cell
+                                     * and the given finite element.
+                                     *
+                                     * The constructor needs a boundary object
+                                     * passed, since this class needs to know
+                                     * how to handle faces which are located
+                                     * on the boundary of the domain. In that
+                                     * case, faces may be curved and the
+                                     * calculation of quadrature points,
+                                     * gradients and the like may need
+                                     * additional effort, depending on the
+                                     * mapping from the unit to the real cell
+                                     * (linear mappings use straight boundary
+                                     * segments, but higher order elements
+                                     * may use other ways.)
+                                     */
+    void reinit (const Triangulation<dim>::face_iterator &,
+                const FiniteElement<dim> &,
+                const Boundary<dim> &);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The index " << arg1
+                   << " is out of range, it should be less than " << arg2);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcAccessToUninitializedField);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcCannotInitializeField);
+    
+  private:
+                                    /**
+                                     * Store the values of the shape functions
+                                     * at the quadrature points. Rows in this
+                                     * matrix denote the values of a single
+                                     * shape function at the different points,
+                                     * columns are for a single point with the
+                                     * different shape functions.
+                                     */
+    dFMatrix             shape_values;
+
+                                    /**
+                                     * Store the gradients of the shape
+                                     * functions at the quadrature points.
+                                     * Since unfortunately the full matrix
+                                     * classes of DEAL are not templated,
+                                     * we have to store them in an
+                                     * archetypic style.
+                                     *
+                                     * This field is reset each time
+                                     * #reinit# is called and contains the
+                                     * gradients on the real element, rather
+                                     * than on the reference element.
+                                     */
+    vector<vector<Point<dim> > >  shape_gradients;
+
+                                    /**
+                                     * Store the gradients of the shape
+                                     * functions at the quadrature points on
+                                     * the unit cell.
+                                     * This field is set up upon construction
+                                     * of the object and contains the gradients
+                                     * on the reference element.
+                                     */
+    vector<vector<Point<dim> > >   unit_shape_gradients;
+    
+                                    /**
+                                     * Store an array of the weights of the
+                                     * quadrature points. This array is
+                                     * set up upon construction.
+                                     */
+    vector<double>       weights;
+
+                                    /**
+                                     * Store an array of weights times the
+                                     * Jacobi determinant at the quadrature
+                                     * points. This function is reset each time
+                                     * #reinit# is called. The Jacobi determinant
+                                     * is actually the reciprocal value of the
+                                     * Jacobi matrices stored in this class,
+                                     * see the general documentation of this
+                                     * class for more information.
+                                     */
+    vector<double>       JxW_values;
+
+                                    /**
+                                     * Array of quadrature points. This array
+                                     * is set up upon calling #reinit# and
+                                     * contains the quadrature points on the
+                                     * real element, rather than on the
+                                     * reference element.
+                                     */
+    vector<Point<dim> >  quadrature_points;
+
+                                    /**
+                                     * Array of quadrature points in the unit
+                                     * cell. This array is set up upon
+                                     * construction and contains the quadrature
+                                     * points on the reference element.
+                                     */
+    vector<Point<dim-1> > unit_quadrature_points;
+    
+                                    /**
+                                     * Array of points denoting the off-point
+                                     * of the ansatz functions. In real space
+                                     * (no-one seems to need the off-point
+                                     * on the unit cell, so no function is
+                                     * provided for this).
+                                     */
+    vector<Point<dim> >  ansatz_points;
+    
+                                    /**
+                                     * Store the jacobi matrices at the
+                                     * different quadrature points. This field
+                                     * is set each time #reinit# is called.
+                                     */
+    vector<dFMatrix>     jacobi_matrices;
+
+                                    /**
+                                     * Store a pointer to the object describing
+                                     * the boundary of the domain.
+                                     */
+    const Boundary<dim> &boundary;
+
+                                    /**
+                                     * Store which fields are to be updated by
+                                     * the reinit function.
+                                     */
+    UpdateFields         update_flags;
+};
+
+
+
+
 /**
   Base class for finite elements in arbitrary dimensions.
   */
@@ -490,7 +812,8 @@ class FiniteElementBase {
                                 vector<Point<dim> > &ansatz_points,
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
-                                const bool           compute_q_points) const;
+                                const bool           compute_q_points,
+                                const Boundary<dim> &boundary) const;
 
                                     /**
                                      * Return the ansatz points this FE has
@@ -498,13 +821,15 @@ class FiniteElementBase {
                                      * given face as a side. This function is
                                      * needed for higher order elements, if
                                      * we want to use curved boundary
-                                     * approximations.
+                                     * 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;
     
                                     /**
@@ -719,7 +1044,8 @@ class FiniteElement<1> : public FiniteElementBase<1> {
                                 vector<Point<1> > &ansatz_points,
                                 const bool         compute_ansatz_points,
                                 vector<Point<1> > &q_points,
-                                const bool         compute_q_points) const;
+                                const bool         compute_q_points,
+                                const Boundary<1> &boundary) const;
 
                                     /**
                                      * Return the ansatz points this FE has
@@ -733,6 +1059,7 @@ class FiniteElement<1> : public FiniteElementBase<1> {
                                      * 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;
 };
 
@@ -869,7 +1196,8 @@ class FiniteElement<2> : public FiniteElementBase<2> {
                                 vector<Point<2> > &ansatz_points,
                                 const bool         compute_ansatz_points,
                                 vector<Point<2> > &q_points,
-                                const bool         compute_q_points) const;
+                                const bool         compute_q_points,
+                                const Boundary<2> &boundary) const;
 
                                     /**
                                      * Return the ansatz points this FE has
@@ -884,6 +1212,7 @@ class FiniteElement<2> : public FiniteElementBase<2> {
                                      * elements.
                                      */
     virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face,
+                                    const Boundary<2>  &boundary,
                                     vector<Point<2> >  &ansatz_points) const;
 };
 
index e4da5ae11cffe16c365c222eac83a11a209f2544..792afc5073e496c3a27d52dc443b6ba3e5cae835 100644 (file)
@@ -80,7 +80,8 @@ class FELinear : public FiniteElement<dim> {
                                 vector<Point<dim> > &ansatz_points,
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
-                                const bool           compute_q_points) const;
+                                const bool           compute_q_points,
+                                const Boundary<dim> &boundary) const;
 
                                     /**
                                      * Return the ansatz points this FE has
@@ -95,6 +96,7 @@ class FELinear : public FiniteElement<dim> {
                                      * elements.
                                      */
     virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+                                    const Boundary<dim>  &boundary,
                                     vector<Point<dim> >  &ansatz_points) const;
 };
 
@@ -164,7 +166,8 @@ class FEQuadratic : public FiniteElement<dim> {
                                 vector<Point<dim> > &ansatz_points,
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
-                                const bool           compute_q_points) const;
+                                const bool           compute_q_points,
+                                const Boundary<dim> &boundary) const;
 };
 
 
@@ -233,7 +236,8 @@ class FECubic : public FiniteElement<dim> {
                                 vector<Point<dim> > &ansatz_points,
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
-                                const bool           compute_q_points) const;
+                                const bool           compute_q_points,
+                                const Boundary<dim> &boundary) const;
 };
 
 
index f8c5f7888809805dcd7a8361252a4569568665d4..fa1ad5121f4102a9957a8347ffa49fc508baa8b7 100644 (file)
@@ -106,7 +106,8 @@ struct AssemblerData {
                   dVector                  &rhs_vector,
                   const Quadrature<dim>    &quadrature,
                   const FiniteElement<dim> &fe,
-                  const UpdateFields       &update_flags);
+                  const UpdateFields       &update_flags,
+                  const Boundary<dim>      &boundary);
     
                                     /**
                                      * Pointer to the dof handler object
@@ -151,7 +152,16 @@ struct AssemblerData {
                                      * FEValues object need to be reinitialized
                                      * on each cell.
                                      */
-    const UpdateFields update_flags;
+    const UpdateFields       update_flags;
+
+                                    /**
+                                     * Store a pointer to the object describing
+                                     * the boundary of the domain. This is
+                                     * necessary, since we may want to use
+                                     * curved faces of cells at the boundary
+                                     * when using higher order elements.
+                                     */
+    const Boundary<dim>     &boundary;
 };
 
 
@@ -229,13 +239,13 @@ class Assembler : public DoFCellAccessor<dim> {
                                      * Pointer to the matrix to be assembled
                                      * by this object.
                                      */
-    dSMatrix               &matrix;
+    dSMatrix         &matrix;
 
                                     /**
                                      * Pointer to the vector to be assembled
                                      * by this object.
                                      */
-    dVector                &rhs_vector;
+    dVector          &rhs_vector;
 
                                     /**
                                      * Pointer to the finite element used for
@@ -258,6 +268,15 @@ class Assembler : public DoFCellAccessor<dim> {
                                      * quadrature points.
                                      */
     FEValues<dim>     fe_values;
+
+                                    /**
+                                     * Store a pointer to the object describing
+                                     * the boundary of the domain. This is
+                                     * necessary, since we may want to use
+                                     * curved faces of cells at the boundary
+                                     * when using higher order elements.
+                                     */
+    const Boundary<dim> &boundary;
 };
 
     
index 0b943cfd888ce9b9991da1e176056bbe0480604a..b4cf45b0365adfec10c063009e28644ff01f083c 100644 (file)
@@ -20,7 +20,8 @@ template <int dim> class DataOut;
 template <int dim> class Function;
 template <int dim> class Equation;
 template <int dim> class Assembler;
-
+template <int dim> class Boundary;
+template <int dim> class StraightBoundary;
 
 
 
@@ -78,6 +79,13 @@ enum NormType {
       induced by hanging nodes.
   \end{itemize}
 
+  The #assemble# function needs an object describing the boundary of the domain,
+  since for higher order finite elements, we may be tempted to use curved faces
+  of cells for better approximation of the boundary. In this case, the
+  transformation from the unit cell to the real cell requires knowledge of
+  the exact boundary of the domain. By default it is assumed that the
+  boundary be approximated by straight segments.
+  
 
   {\bf Solving}
 
@@ -294,7 +302,8 @@ class ProblemBase {
                           const Quadrature<dim>    &q,
                           const FiniteElement<dim> &fe,
                           const UpdateFields       &update_flags,
-                          const DirichletBC        &dirichlet_bc = DirichletBC());
+                          const DirichletBC        &dirichlet_bc = DirichletBC(),
+                          const Boundary<dim>      &boundary = StraightBoundary<dim>());
     
                                     /**
                                      * Solve the system of equations.
@@ -315,7 +324,8 @@ class ProblemBase {
                               dVector                  &difference,
                               const Quadrature<dim>    &q,
                               const FiniteElement<dim> &fe,
-                              const NormType           &norm) const;
+                              const NormType           &norm,
+                              const Boundary<dim> &boundary=StraightBoundary<dim>()) const;
     
                                     /**
                                      * Initialize the #DataOut# object with
@@ -358,6 +368,7 @@ class ProblemBase {
                                      */
     virtual void make_boundary_value_list (const DirichletBC        &dirichlet_bc,
                                           const FiniteElement<dim> &fe,
+                                          const Boundary<dim>      &boundary,
                                           map<int,double>          &boundary_values) const;
     
                                     /**
@@ -433,7 +444,8 @@ class ProblemBase {
                             dVector           &solution,
                             dVector           &right_hand_side,
                             const DirichletBC &dirichlet_bc,
-                            const FiniteElement<dim> &fe);
+                            const FiniteElement<dim> &fe,
+                            const Boundary<dim>      &boundary);
     
     friend class Assembler<dim>;
 };
index b28f84258421253a0f00d19773ad0e0a40002f80..911f77d4385874dc4eb76ede989b2410a6f394e8 100644 (file)
@@ -4,7 +4,7 @@
 #include <fe/quadrature.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_accessor.h>
-
+#include <grid/tria_boundary.h>
 
 
 
@@ -52,8 +52,8 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
 template <int dim>
 double FEValues<dim>::shape_value (const unsigned int i,
                                   const unsigned int j) const {
-  Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
-  Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+  Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+  Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
 
   return shape_values(i,j);
 };
@@ -64,8 +64,8 @@ template <int dim>
 const Point<dim> &
 FEValues<dim>::shape_grad (const unsigned int i,
                           const unsigned int j) const {
-  Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
-  Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+  Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+  Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
   Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
 
   return shape_gradients[i][j];
@@ -105,7 +105,8 @@ 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 FiniteElement<dim>                         &fe,
+                           const Boundary<dim>                              &boundary) {
                                   // fill jacobi matrices and real
                                   // quadrature points
   if ((update_flags | update_jacobians) ||
@@ -117,7 +118,8 @@ void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &ce
                       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
@@ -149,8 +151,7 @@ void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &ce
                                   // determinant
   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();
     };
@@ -239,7 +240,8 @@ void FiniteElementBase<dim>::fill_fe_values (const typename Triangulation<dim>::
                                             vector<Point<dim> > &,
                                             const bool,
                                             vector<Point<dim> > &,
-                                            const bool) const {
+                                            const bool,
+                                            const Boundary<dim> &) const {
   Assert (false, ExcPureFunctionCalled());
 };
 
@@ -247,6 +249,7 @@ void FiniteElementBase<dim>::fill_fe_values (const typename Triangulation<dim>::
 
 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());
 };
@@ -271,7 +274,8 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce
                                       vector<Point<1> > &ansatz_points,
                                       const bool         compute_ansatz_points,
                                       vector<Point<1> > &q_points,
-                                      const bool         compute_q_points) const {
+                                      const bool         compute_q_points,
+                                      const Boundary<1> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -313,6 +317,7 @@ 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?
   Assert (false, ExcPureFunctionCalled());
@@ -336,13 +341,15 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                       vector<Point<2> > &,
                                       const bool,
                                       vector<Point<2> > &,
-                                      const bool) const {
+                                      const bool,
+                                      const Boundary<2> &) const {
   Assert (false, ExcPureFunctionCalled());
 };
 
 
 
 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());
index 0604e9e7a2a9589ed43e1fe60069f9e0999de65d..5eafa7b819963012ed39ba43f64b35cbcf9a6bdd 100644 (file)
@@ -85,17 +85,19 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                  vector<Point<1> > &ansatz_points,
                                  const bool         compute_ansatz_points,
                                  vector<Point<1> > &q_points,
-                                 const bool         compute_q_points) const {
+                                 const bool         compute_q_points,
+                                 const Boundary<1> &boundary) const {
                                   // simply pass down
   FiniteElement<1>::fill_fe_values (cell, unit_points,
                                    jacobians, compute_jacobians,
                                    ansatz_points, compute_ansatz_points,
-                                   q_points, compute_q_points);
+                                   q_points, compute_q_points, boundary);
 };
 
 
 
 void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &,
+                                     const Boundary<1> &,
                                      vector<Point<1> > &) const {
   Assert (false, ExcPureFunctionCalled());
 };
@@ -246,7 +248,8 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                  vector<Point<2> > &ansatz_points,
                                  const bool         compute_ansatz_points,
                                  vector<Point<2> > &q_points,
-                                 const bool         compute_q_points) const {
+                                 const bool         compute_q_points,
+                                 const Boundary<2> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -335,6 +338,7 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
 
 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 {
   Assert (ansatz_points.size() == (1<<(dim-1)),
          typename FiniteElementBase<dim>::ExcWrongFieldDimension (ansatz_points.size(),
@@ -362,12 +366,13 @@ void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell
                                     vector<Point<1> > &ansatz_points,
                                     const bool         compute_ansatz_points,
                                     vector<Point<1> > &q_points,
-                                    const bool         compute_q_points) const {
+                                    const bool         compute_q_points,
+                                    const Boundary<1> &boundary) const {
                                   // simply pass down
   FiniteElement<1>::fill_fe_values (cell, unit_points,
                                    jacobians, compute_jacobians,
                                    ansatz_points, compute_ansatz_points,
-                                   q_points, compute_q_points);
+                                   q_points, compute_q_points, boundary);
 };
 
 
@@ -418,7 +423,8 @@ void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                     vector<Point<2> > &ansatz_points,
                                     const bool,
                                     vector<Point<2> > &q_points,
-                                    const bool) const {
+                                    const bool,
+                                    const Boundary<2> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
@@ -446,12 +452,13 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                 vector<Point<1> > &ansatz_points,
                                 const bool         compute_ansatz_points,
                                 vector<Point<1> > &q_points,
-                                const bool         compute_q_points) const {
+                                const bool         compute_q_points,
+                                const Boundary<1> &boundary) const {
                                   // simply pass down
   FiniteElement<1>::fill_fe_values (cell, unit_points,
                                    jacobians, compute_jacobians,
                                    ansatz_points, compute_ansatz_points,
-                                   q_points, compute_q_points);
+                                   q_points, compute_q_points, boundary);
 };
 
 
@@ -494,7 +501,8 @@ void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                 vector<Point<2> > &ansatz_points,
                                 const bool,
                                 vector<Point<2> > &q_points,
-                                const bool) const {
+                                const bool,
+                                const Boundary<2> &) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
   Assert (q_points.size() == unit_points.size(),
index cf47c72386572bd60e1b1b067d2de987c52a592c..616ac3b0abcdf311b1b044b6902044f46c3f835e 100644 (file)
@@ -51,7 +51,8 @@ AssemblerData<dim>::AssemblerData (const DoFHandler<dim>    &dof,
                                   dVector                  &rhs_vector,
                                   const Quadrature<dim>    &quadrature,
                                   const FiniteElement<dim> &fe,
-                                  const UpdateFields       &update_flags) :
+                                  const UpdateFields       &update_flags,
+                                  const Boundary<dim>      &boundary) :
                dof(dof),
                assemble_matrix(assemble_matrix),
                assemble_rhs(assemble_rhs),
@@ -59,7 +60,8 @@ AssemblerData<dim>::AssemblerData (const DoFHandler<dim>    &dof,
                rhs_vector(rhs_vector),
                quadrature(quadrature),
                fe(fe),
-               update_flags(update_flags) {};
+               update_flags(update_flags),
+               boundary(boundary) {};
 
 
 
@@ -80,7 +82,8 @@ Assembler<dim>::Assembler (Triangulation<dim> *tria,
                fe(((AssemblerData<dim>*)local_data)->fe),
                fe_values (((AssemblerData<dim>*)local_data)->fe,
                           ((AssemblerData<dim>*)local_data)->quadrature,
-                          ((AssemblerData<dim>*)local_data)->update_flags)
+                          ((AssemblerData<dim>*)local_data)->update_flags),
+               boundary(((AssemblerData<dim>*)local_data)->boundary)
 {
   Assert (matrix.m() == dof_handler->n_dofs(), ExcInvalidData());
   Assert (matrix.n() == dof_handler->n_dofs(), ExcInvalidData());
@@ -98,7 +101,8 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
   fe_values.reinit (Triangulation<dim>::cell_iterator (tria,
                                                       present_level,
                                                       present_index),
-                   fe);
+                   fe,
+                   boundary);
   const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs;
 
                                   // clear cell matrix
index e5b3cd1875f7d1a551dbc47fadd134ea86fc6f62..6e0aff2d484d2f41768b1761881b6765eeedce43 100644 (file)
@@ -77,7 +77,8 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                                 const Quadrature<dim>    &quadrature,
                                 const FiniteElement<dim> &fe,
                                 const UpdateFields       &update_flags,
-                                const DirichletBC        &dirichlet_bc) {
+                                const DirichletBC        &dirichlet_bc,
+                                const Boundary<dim>      &boundary) {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
   
   system_sparsity.reinit (dof_handler->n_dofs(),
@@ -105,7 +106,8 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                           right_hand_side,
                           quadrature,
                           fe,
-                          update_flags);
+                          update_flags,
+                          boundary);
   active_assemble_iterator assembler (tria,
                                      tria->begin_active()->level(),
                                      tria->begin_active()->index(),
@@ -127,7 +129,7 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                                   // in the docs
   apply_dirichlet_bc (system_matrix, solution,
                      right_hand_side,
-                     dirichlet_bc, fe);
+                     dirichlet_bc, fe, boundary);
   
 };
 
@@ -158,7 +160,8 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
                                             dVector                  &difference,
                                             const Quadrature<dim>    &q,
                                             const FiniteElement<dim> &fe,
-                                            const NormType           &norm) const {
+                                            const NormType           &norm,
+                                            const Boundary<dim>      &boundary) const {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());  
   Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
 
@@ -187,7 +190,7 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
     {
       double diff=0;
                                       // initialize for this cell
-      fe_values.reinit (tria_cell, fe);
+      fe_values.reinit (tria_cell, fe, boundary);
 
       switch (norm) 
        {
@@ -375,7 +378,8 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
                                           dVector  &solution,
                                           dVector  &right_hand_side,
                                           const DirichletBC &dirichlet_bc,
-                                          const FiniteElement<dim> &fe) {
+                                          const FiniteElement<dim> &fe,
+                                          const Boundary<dim>      &boundary) {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
   Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
          ExcInvalidBoundaryIndicator());
@@ -387,7 +391,7 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
                                   // the lines in 2D being subject to
                                   // different bc's), the last value is taken
   map<int,double> boundary_values;
-  make_boundary_value_list (dirichlet_bc, fe, boundary_values);
+  make_boundary_value_list (dirichlet_bc, fe, boundary, boundary_values);
 
   map<int,double>::const_iterator dof, endd;
   const unsigned int n_dofs   = (unsigned int)matrix.m();
@@ -477,6 +481,7 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
 void
 ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
                                          const FiniteElement<1> &,
+                                         const Boundary<1> &,
                                          map<int,double>   &) const {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());  
   Assert (false, ExcNotImplemented());
@@ -487,8 +492,9 @@ ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
 
 template <int dim>
 void
-ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
+ProblemBase<dim>::make_boundary_value_list (const DirichletBC        &dirichlet_bc,
                                            const FiniteElement<dim> &fe,
+                                           const Boundary<dim>      &boundary,
                                            map<int,double>   &boundary_values) const {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
 
@@ -521,7 +527,7 @@ ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
        face_dofs.erase (face_dofs.begin(), face_dofs.end());
        dof_values.erase (dof_values.begin(), dof_values.end());
        face->get_dof_indices (face_dofs);
-       fe.face_ansatz_points (tface, dof_locations);
+       fe.face_ansatz_points (tface, boundary, dof_locations);
        (*function_ptr).second->value_list (dof_locations, dof_values);
 
                                         // enter into list

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.