]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Full support for Dirichlet BC and bug fixes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 3 Apr 1998 17:44:30 +0000 (17:44 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 3 Apr 1998 17:44:30 +0000 (17:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@134 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.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/base.cc

index 20aebc1c9f93d2b1d9ecf589a6364a4584c435d4..1d39fb5de81bd7086a1f3c8bceeb4c564d7677aa 100644 (file)
@@ -1095,7 +1095,14 @@ class DoFHandler : public DoFDimensionInfo<dim> {
 
                                     /**
                                      * Store a copy of the finite element given
-                                     * latest for the distribution of dofs.
+                                     * latest for the distribution of dofs. In
+                                     * fact, since the FE class itself has not
+                                     * much functionality, this object only
+                                     * stores numbers such as the number of
+                                     * dofs per line etc. Calling any of the
+                                     * more specific functions will result in
+                                     * an error about calling a pure virtual
+                                     * function.
                                      */
     FiniteElement<dim>       *selected_fe;
 
index 592dbc5ba3185baa3e085db67e31ac0ff00ac436..f1188b0b49844cb6813e131cc6394e8f98830e30 100644 (file)
@@ -35,7 +35,7 @@ template <int dim> class Quadrature;
   the real cell at the quadrature points and so on.
 
   The Jacobian matrix is defined to be
-  $$ J_{ij} = {d\xi_i \over d\x_j} $$
+  $$ J_{ij} = {d\xi_i \over dx_j} $$
   which is the form needed to compute the gradient on the real cell from
   the gradient on the unit cell. If we want to transform the area element
   $dx dy$ from the real to the unit cell, we have to take the determinant of
@@ -392,7 +392,7 @@ class FiniteElementBase {
                                      * 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 ommit the
+                                     * this function. Therefore, we omit the
                                      * #=0# signature and implement this function
                                      * by throwing an exception.
                                      * #p# is a point on the reference element.
@@ -406,7 +406,7 @@ class FiniteElementBase {
                                      * 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 ommit the
+                                     * this function. Therefore, we omit the
                                      * #=0# signature and implement this function
                                      * by throwing an exception.
                                      * #p# is a point on the reference element,
@@ -472,6 +472,10 @@ class FiniteElementBase {
                                      * 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,
@@ -481,6 +485,21 @@ class FiniteElementBase {
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points) 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.
+                                     *
+                                     * The function assumes that the fields
+                                     * already have the right number of
+                                     * elements.
+                                     */
+    virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+                                    vector<Point<dim> >  &ansatz_points) const;
     
                                     /**
                                      * Comparison operator. We also check for
@@ -512,7 +531,13 @@ class FiniteElementBase {
                                      * Exception
                                      */
     DeclException0 (ExcNotImplemented);
-    
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcWrongFieldDimension,
+                   int, int,
+                   << "The field has not the assumed dimension " << arg2
+                   << ", but has " << arg1 << " elements.");
   protected:
                                     /**
                                      * Have #N=2^dim# matrices keeping the
@@ -598,6 +623,12 @@ class FiniteElement;
   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:
@@ -670,6 +701,10 @@ class FiniteElement<1> : public FiniteElementBase<1> {
                                      * 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,
@@ -679,6 +714,20 @@ class FiniteElement<1> : public FiniteElementBase<1> {
                                 const bool         compute_ansatz_points,
                                 vector<Point<1> > &q_points,
                                 const bool         compute_q_points) 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,
+                                    vector<Point<1> >  &ansatz_points) const;
 };
 
 
@@ -719,6 +768,12 @@ class FiniteElement<1> : public FiniteElementBase<1> {
 
   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> {
   public:
@@ -796,6 +851,10 @@ class FiniteElement<2> : public FiniteElementBase<2> {
                                      * 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<2>::cell_iterator &cell,
                                 const vector<Point<2> >               &unit_points,
@@ -805,6 +864,21 @@ class FiniteElement<2> : public FiniteElementBase<2> {
                                 const bool         compute_ansatz_points,
                                 vector<Point<2> > &q_points,
                                 const bool         compute_q_points) 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.
+                                     *
+                                     * The function assumes that the fields
+                                     * already have the right number of
+                                     * elements.
+                                     */
+    virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face,
+                                    vector<Point<2> >  &ansatz_points) const;
 };
 
 
index c6171f526864fc3af4b9932996a135e55910ab48..e4da5ae11cffe16c365c222eac83a11a209f2544 100644 (file)
@@ -68,6 +68,10 @@ class FELinear : public FiniteElement<dim> {
                                      * 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,
@@ -77,6 +81,21 @@ class FELinear : public FiniteElement<dim> {
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points) 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.
+                                     */
+    virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+                                    vector<Point<dim> >  &ansatz_points) const;
 };
 
 
@@ -133,6 +152,10 @@ class FEQuadratic : public FiniteElement<dim> {
                                      * 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,
@@ -198,6 +221,10 @@ class FECubic : public FiniteElement<dim> {
                                      * 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,
index 7f0422227a09c74816cda292650b1e4d965eecbd..4a716bb56deaf27bf856fd6bb68efe018fe4e2f4 100644 (file)
@@ -331,8 +331,9 @@ class ProblemBase {
                                      * See the general doc for more
                                      * information.
                                      */
-    virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
-                                          map<int,double>   &boundary_values) const;
+    virtual void make_boundary_value_list (const DirichletBC        &dirichlet_bc,
+                                          const FiniteElement<dim> &fe,
+                                          map<int,double>          &boundary_values) const;
     
                                     /**
                                      * Exception
@@ -406,7 +407,8 @@ class ProblemBase {
     void apply_dirichlet_bc (dSMatrix          &matrix,
                             dVector           &solution,
                             dVector           &right_hand_side,
-                            const DirichletBC &dirichlet_bc);
+                            const DirichletBC &dirichlet_bc,
+                            const FiniteElement<dim> &fe);
     
     friend class Assembler<dim>;
 };
index 8543071a75e970d89f6885ade9bde79fda887fc5..c4907af43f7e091a890b413a0e4615a1c18ac1b8 100644 (file)
@@ -37,6 +37,7 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                JxW_values(quadrature.n_quadrature_points, 0),
                quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
                unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
+               ansatz_points (fe.total_dofs, Point<dim>()),
                jacobi_matrices (quadrature.n_quadrature_points,
                                 dFMatrix(dim,dim)),
                update_flags (update_flags)
@@ -253,6 +254,14 @@ 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 &,
+                                                vector<Point<dim> > &) const {
+  Assert (false, ExcPureFunctionCalled());
+};
+  
+
+
 /*------------------------------- FiniteElement ----------------------*/
 
 
@@ -272,6 +281,14 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce
                                       const bool         compute_ansatz_points,
                                       vector<Point<1> > &q_points,
                                       const bool         compute_q_points) 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 (ansatz_points.size() == total_dofs,
+         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
+
                                   // local mesh width
   const double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
 
@@ -306,6 +323,13 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce
 
 
 
+void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &,
+                                          vector<Point<1> > &) const {
+                                  // is this function useful in 1D?
+  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) &&
@@ -328,6 +352,14 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
 
 
 
+void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &,
+                                          vector<Point<2> > &) const {
+                                  // is this function useful in 1D?
+  Assert (false, ExcPureFunctionCalled());
+};
+
+
+
 
 
 /*------------------------------- Explicit Instantiations -------------*/
index 93cb12965e89915a5bfc17d22d7b1bb50ca3327f..0604e9e7a2a9589ed43e1fe60069f9e0999de65d 100644 (file)
@@ -5,6 +5,8 @@
 #include <grid/tria_accessor.h>
 
 
+
+
 FELinear<1>::FELinear () :
                FiniteElement<1> (1, 0)
 {
@@ -93,6 +95,12 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
 
 
 
+void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &,
+                                     vector<Point<1> > &) const {
+  Assert (false, ExcPureFunctionCalled());
+};
+
+
 
 
 FELinear<2>::FELinear () :
@@ -239,8 +247,15 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                  const bool         compute_ansatz_points,
                                  vector<Point<2> > &q_points,
                                  const bool         compute_q_points) 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 (ansatz_points.size() == total_dofs,
+         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+  
   const unsigned int dim=2;
-  const unsigned int n_vertices=4;
+  const unsigned int n_vertices=(1<<dim);
   unsigned int n_points=unit_points.size();
 
   Point<dim> vertices[n_vertices];
@@ -318,6 +333,20 @@ 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,
+                                       vector<Point<dim> >  &ansatz_points) const {
+  Assert (ansatz_points.size() == (1<<(dim-1)),
+         typename FiniteElementBase<dim>::ExcWrongFieldDimension (ansatz_points.size(),
+                                                                  1<<(dim-1)));
+  
+  for (unsigned int vertex=0; vertex<(1<<(dim-1)); ++vertex)
+    ansatz_points[vertex] = face->vertex(vertex);
+};
+
+
+
+
 
 
 
@@ -383,13 +412,20 @@ FEQuadratic<dim>::shape_grad (const unsigned int i,
 
 
 void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
-                                    const vector<Point<2> >               &,
-                                    vector<dFMatrix>  &,
+                                    const vector<Point<2> >               &unit_points,
+                                    vector<dFMatrix>  &jacobians,
                                     const bool,
-                                    vector<Point<2> > &,
+                                    vector<Point<2> > &ansatz_points,
                                     const bool,
-                                    vector<Point<2> > &,
+                                    vector<Point<2> > &q_points,
                                     const bool) 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 (ansatz_points.size() == total_dofs,
+         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
   Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
 };
 
@@ -452,13 +488,20 @@ FECubic<dim>::shape_grad (const unsigned int i,
 
 
 void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
-                                const vector<Point<2> >               &,
-                                vector<dFMatrix>  &,
+                                const vector<Point<2> >               &unit_points,
+                                vector<dFMatrix>  &jacobians,
                                 const bool,
-                                vector<Point<2> > &,
+                                vector<Point<2> > &ansatz_points,
                                 const bool,
-                                vector<Point<2> > &,
+                                vector<Point<2> > &q_points,
                                 const bool) 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 (ansatz_points.size() == total_dofs,
+         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
   Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
 };
 
index c3c3edbb9412db4f785dba08205c37c64b5b9d4e..7b555ba9164d5832294d80b7293a06f18dfad706 100644 (file)
@@ -113,7 +113,8 @@ void ProblemBase<dim>::assemble (const Equation<dim>               &equation,
                                   // apply Dirichlet bc as described
                                   // in the docs
   apply_dirichlet_bc (system_matrix, solution,
-                     right_hand_side, dirichlet_bc);
+                     right_hand_side,
+                     dirichlet_bc, fe);
   
                                   // condense system matrix in-place
   constraints.condense (system_matrix);
@@ -307,11 +308,12 @@ template <int dim>
 void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
                                           dVector  &solution,
                                           dVector  &right_hand_side,
-                                          const DirichletBC &dirichlet_bc) {
+                                          const DirichletBC &dirichlet_bc,
+                                          const FiniteElement<dim> &fe) {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
   Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
          ExcInvalidBoundaryIndicator());
-  
+
                                   // first make up a list of dofs subject
                                   // to any boundary condition and which
                                   // value they take; if a node occurs
@@ -319,7 +321,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, boundary_values);
+  make_boundary_value_list (dirichlet_bc, fe, boundary_values);
 
   map<int,double>::const_iterator dof, endd;
   const unsigned int n_dofs   = (unsigned int)matrix.m();
@@ -346,10 +348,16 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
                                       // zero, take the first main
                                       // diagonal entry we can find, or
                                       // one if no nonzero main diagonal
-                                      // element exists.
+                                      // element exists. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      double new_rhs;
       if (matrix.diag_element((*dof).first) != 0.0)
-       right_hand_side((*dof).first) = (*dof).second *
-                                       matrix.diag_element((*dof).first);
+       new_rhs = right_hand_side((*dof).first)
+               = (*dof).second * matrix.diag_element((*dof).first);
       else
        {
          double first_diagonal_entry = 1;
@@ -362,19 +370,19 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
          
          matrix.set((*dof).first, (*dof).first,
                     first_diagonal_entry);
-         right_hand_side((*dof).first) = (*dof).second * first_diagonal_entry;
+         new_rhs = right_hand_side((*dof).first)
+                 = (*dof).second * first_diagonal_entry;
        };
-
+      
                                       // store the only nonzero entry
                                       // of this line for the Gauss
                                       // elimination step
       const double diagonal_entry = matrix.diag_element((*dof).first);
-      
 
                                       // do the Gauss step
-      for (unsigned int row=0; row<n_dofs; ++row)
+      for (unsigned int row=0; row<n_dofs; ++row) 
        for (int j=sparsity_rowstart[row];
-            j<sparsity_rowstart[row+1]; ++row)
+            j<sparsity_rowstart[row+1]; ++j)
          if ((sparsity_colnums[j] == (signed int)(*dof).first) &&
              ((signed int)row != (*dof).first))
                                             // this line has an entry
@@ -384,7 +392,7 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
            {
                                               // correct right hand side
              right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
-                                     (*dof).second;
+                                     new_rhs;
              
                                               // set matrix entry to zero
              matrix.global_entry(j) = 0.;
@@ -402,6 +410,7 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
 
 void
 ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
+                                         const FiniteElement<1> &,
                                          map<int,double>   &) const {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());  
   Assert (false, ExcNotImplemented());
@@ -413,13 +422,28 @@ ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
 template <int dim>
 void
 ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
+                                           const FiniteElement<dim> &fe,
                                            map<int,double>   &boundary_values) const {
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
-  
+
+                                  // use two face iterators, since we need
+                                  // a DoF-iterator for the dof indices, but
+                                  // a Tria-iterator for the fe object
   DoFHandler<dim>::active_face_iterator face = dof_handler->begin_active_face(),
                                        endf = dof_handler->end_face();
+  Triangulation<dim>::active_face_iterator tface = tria->begin_active_face();
+  
   DirichletBC::const_iterator function_ptr;
-  for (; face!=endf; ++face)
+
+                                  // field to store the indices of dofs
+                                  // initialize once to get the size right
+                                  // for the following fields.
+  vector<int>         face_dofs;
+  face->get_dof_indices (face_dofs);
+  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+  vector<double>      dof_values;
+       
+  for (; face!=endf; ++face, ++tface)
     if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
        dirichlet_bc.end()) 
                                       // face is subject to one of the
@@ -428,17 +452,10 @@ ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
                                         // get indices, physical location and
                                         // boundary values of dofs on this
                                         // face
-       vector<int>         face_dofs;
-       vector<Point<dim> > dof_locations;
-       vector<double>      dof_values;
-       
+       face_dofs.erase (face_dofs.begin(), face_dofs.end());
+       dof_values.erase (dof_values.begin(), dof_values.end());
        face->get_dof_indices (face_dofs);
-
-//physical location
-//     Assert (false, ExcNotImplemented());
-       dof_locations.insert (dof_locations.begin(),
-                             face_dofs.size(), Point<dim>());
-       
+       fe.face_ansatz_points (tface, 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.