]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added a second constructor to FE_Q that allows to generate it based on support points...
authorkormann <kormann@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 7 Sep 2008 19:44:38 +0000 (19:44 +0000)
committerkormann <kormann@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 7 Sep 2008 19:44:38 +0000 (19:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@16756 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/quadrature_lib.cc
deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/doc/authors.html
deal.II/doc/news/changes.h

index 05d3d2ec3654663a94925e075d1c55b062b74a65..72fd24800875876495c7cd9fb0ccdaf7b749d9b0 100644 (file)
@@ -50,6 +50,19 @@ QGauss<0>::QGauss (const unsigned int)
 
 
 
+template <>
+QGaussLobatto<0>::QGaussLobatto (const unsigned int)
+                :
+                Quadrature<0> (0)
+{
+                                   // this function has to be provided to
+                                   // avoid certain linker failures, but it
+                                   // should never be called
+  Assert (false, ExcInternalError());
+}
+
+
+
 template <>
 QGauss<1>::QGauss (const unsigned int n)
                 :
index 0975f0dbed7f25ca8a3ab513dce70cc98196d8d9..35c147ff66d353a4561090c4d6c6866d3516795f 100644 (file)
@@ -24,21 +24,23 @@ DEAL_II_NAMESPACE_OPEN
 /*@{*/
 
 /**
- * Implementation of a scalar Lagrange finite element @p Qp that
- * yields the finite element space of continuous, piecewise polynomials
- * of degree @p p in each coordinate direction. This class is realized
- * using tensor product polynomials based on equidistant support
- * points.
+ * Implementation of a scalar Lagrange finite element @p Qp that yields the
+ * finite element space of continuous, piecewise polynomials of degree @p p in
+ * each coordinate direction. This class is realized using tensor product
+ * polynomials based on equidistant or given support points.
  *
- * The constructor of this class takes the degree @p p of this finite
- * element.
+ * The standard constructor of this class takes the degree @p p of this finite
+ * element. Alternatively, it can take a quadrature formula @p points defining
+ * the support points of the Lagrange interpolation in one coordinate direction.
  *
  * <h3>Implementation</h3>
  *
- * The constructor creates a TensorProductPolynomials object
- * that includes the tensor product of @p LagrangeEquidistant
- * polynomials of degree @p p. This @p TensorProductPolynomials
- * object provides all values and derivatives of the shape functions.
+ * The constructor creates a TensorProductPolynomials object that includes the
+ * tensor product of @p LagrangeEquidistant polynomials of degree @p p. This
+ * @p TensorProductPolynomials object provides all values and derivatives of
+ * the shape functions.  In case a quadrature rule is given, the constructure
+ * creates a TensorProductPolynomials object that includes the tensor product
+ * of @p Lagrange polynomials with the support points from @p points.
  *
  * Furthermore the constructor filles the @p interface_constraints,
  * the @p prolongation (embedding) and the @p restriction
@@ -232,7 +234,7 @@ DEAL_II_NAMESPACE_OPEN
  *   @endverbatim
  * </ul>
  *
- * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004
+ * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004; Katharina Kormann, 2008; Martin Kronbichler, 2008
  */
 template <int dim>
 class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
@@ -243,7 +245,19 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      * polynomials of degree @p p.
                                      */
     FE_Q (const unsigned int p);
-    
+       
+                                    /**
+                                     * Constructor for tensor product
+                                     * polynomials with support points @p
+                                     * points based on a one-dimensional
+                                     * quadrature formula. The degree of the
+                                     * finite element is
+                                     * <tt>points.size()-1</tt>.  Note that
+                                     * the first point has to be 0 and the
+                                     * last one 1.
+                                     */
+
+    FE_Q (const Quadrature<1> &points);    
                                     /**
                                      * Return a string that uniquely
                                      * identifies a finite
@@ -494,6 +508,14 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      * from the constructor.
                                      */
     void initialize_constraints ();
+       
+                                            /**
+                                     * Initialize the hanging node
+                                     * constraints matrices. Called from the
+                                     * constructor in case the finite element
+                                     * is based on quadrature points.
+                                     */
+    void initialize_constraints (const Quadrature<1> &points);
 
                                     /**
                                      * Initialize the embedding
@@ -517,6 +539,15 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      * constructor.
                                      */
     void initialize_unit_support_points ();
+       
+                                            /**
+                                     * Initialize the @p unit_support_points
+                                     * field of the FiniteElement
+                                     * class. Called from the constructor in
+                                     * case the finite element is based on
+                                     * quadrature points.
+                                     */
+    void initialize_unit_support_points (const Quadrature<1> &points);
 
                                     /**
                                      * Initialize the
@@ -526,6 +557,15 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      * constructor.
                                      */
     void initialize_unit_face_support_points ();
+       
+                                            /**
+                                     * Initialize the @p
+                                     * unit_face_support_points field of the
+                                     * FiniteElement class. Called from the
+                                     * constructor in case the finite element
+                                     * is based on quadrature points.
+                                     */
+    void initialize_unit_face_support_points (const Quadrature<1> &points);
 
                                     /**
                                      * Initialize the
@@ -571,13 +611,13 @@ std::vector<unsigned int>
 FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
 
 template <>
-void FE_Q<1>::initialize_constraints ();
+void FE_Q<1>::initialize_constraints (const Quadrature<1>&);
 
 template <>
-void FE_Q<2>::initialize_constraints ();
+void FE_Q<2>::initialize_constraints (const Quadrature<1>&);
 
 template <>
-void FE_Q<3>::initialize_constraints ();
+void FE_Q<3>::initialize_constraints (const Quadrature<1>&);
 
 
 DEAL_II_NAMESPACE_CLOSE
index 4b06ab478684b1abc7167f5ff73ce39872f8debc..e0f9bfc95115385f2316fa3873ee67b52855dadb 100644 (file)
@@ -12,6 +12,7 @@
 //---------------------------------------------------------------------------
 
 #include <base/quadrature.h>
+#include <base/quadrature_lib.h>
 #include <base/template_constraints.h>
 #include <fe/fe_dgq.h>
 #include <fe/fe_tools.h>
@@ -231,7 +232,6 @@ FE_DGQ<dim>::get_name () const
 
   std::ostringstream namebuf;  
   namebuf << "FE_DGQ<" << dim << ">(" << this->degree << ")";
-
   return namebuf.str();
 }
 
@@ -666,8 +666,57 @@ FE_DGQArbitraryNodes<dim>::get_name () const
                                   // FE_DGQArbitraryNodes since
                                   // there is no initialization by
                                   // a degree value.
-  std::ostringstream namebuf;  
-  namebuf << "FE_DGQArbitraryNodes<" << dim << ">(" << this->degree << ")";
+  std::ostringstream namebuf;
+
+  bool type = true;  
+  const unsigned int n_points = this->degree +1;
+  std::vector<double> points(n_points);
+  const unsigned int dofs_per_cell = this->dofs_per_cell;
+  const std::vector<Point<dim> > &unit_support_points = this->unit_support_points;
+  unsigned int index = 0;
+
+                                  // Decode the support points
+                                  // in one coordinate direction.
+  for (unsigned int j=0;j<dofs_per_cell;j++)
+    {
+      if ((dim>1) ? (unit_support_points[j](1)==0 && 
+          ((dim>2) ? unit_support_points[j](2)==0: true)) : true)
+       {
+         points[index++] = unit_support_points[j](0);
+       }
+    }
+  Assert (index == n_points,
+         ExcMessage ("Could not decode support points in one coordinate direction."));
+
+                                 // Check whether the support
+                                 // points are equidistant.
+  for(unsigned int j=0;j<n_points;j++)
+    if (std::fabs(points[j] - (double)j/this->degree) > 1e-15)
+      {
+       type = false;
+       break;
+      }
+
+  if (type == true)    
+    namebuf << "FE_DGQ<" << dim << ">(" << this->degree << ")";
+  else
+    {
+
+                                 // Check whether the support
+                                 // points come from QGaussLobatto.
+      const QGaussLobatto<1> points_gl(n_points);
+      type = true;
+      for(unsigned int j=0;j<n_points;j++)
+       if (points[j] != points_gl.point(j)(0))
+         {
+           type = false;
+           break;
+         }
+      if(type == true)
+       namebuf << "FE_DGQArbitraryNodes<" << dim << ">(QGaussLobatto(" << this->degree+1 << "))";
+      else
+       namebuf << "FE_DGQArbitraryNodes<" << dim << ">(QUnknownNodes(" << this->degree << "))";
+    }
 
   return namebuf.str();
 }
index 29aa0d8eca83d093a730df87457262183ec11b40..5377a8fb474045f469d9df285b904c080abf9dd2 100644 (file)
 #include <base/template_constraints.h>
 #include <fe/fe_q.h>
 #include <fe/fe_tools.h>
+#include <base/quadrature_lib.h>
 
 #include <vector>
-
+#include <iostream>
 #include <sstream>
 
 DEAL_II_NAMESPACE_OPEN
@@ -209,6 +210,58 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
 
 
 
+template <int dim>
+FE_Q<dim>::FE_Q (const Quadrature<1> &points)
+               :
+               FE_Poly<TensorProductPolynomials<dim>, dim> (
+                 TensorProductPolynomials<dim>(Polynomials::Lagrange::generate_complete_basis(points.get_points())),
+                 FiniteElementData<dim>(get_dpo_vector(points.n_quadrature_points-1),
+                                        1, points.n_quadrature_points-1,
+                                        FiniteElementData<dim>::H1),
+                 std::vector<bool> (1, false),
+                 std::vector<std::vector<bool> >(1, std::vector<bool>(1,true))),
+               face_index_map(FE_Q_Helper::invert_numbering(face_lexicographic_to_hierarchic_numbering (points.n_quadrature_points-1)))
+{
+  const unsigned int degree = points.n_quadrature_points-1;
+  
+  Assert (degree > 0,
+          ExcMessage ("This element can only be used for polynomial degrees "
+                      "at least zero"));
+  Assert (points.point(0)(0) == 0,
+         ExcMessage ("The first support point has to be zero."));
+  Assert (points.point(degree)(0) == 1,
+         ExcMessage ("The last support point has to be one."));
+  
+  std::vector<unsigned int> renumber (this->dofs_per_cell);
+  FETools::hierarchic_to_lexicographic_numbering (*this, renumber);
+  this->poly_space.set_numbering(renumber);
+  
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points (points);
+  initialize_unit_face_support_points (points);
+
+                                  // compute constraint, embedding
+                                  // and restriction matrices
+  initialize_constraints (points);
+
+                                  // Reinit the vectors of
+                                  // restriction and prolongation
+                                  // matrices to the right sizes
+  this->reinit_restriction_and_prolongation_matrices();
+
+                                  // Fill prolongation matrices with
+                                  // embedding operators
+  FETools::compute_embedding_matrices (*this, this->prolongation);
+
+                                  // Fill restriction matrices with
+                                  // L2-projection
+  FETools::compute_projection_matrices (*this, this->restriction);
+  initialize_quad_dof_index_permutation();
+}
+
+
+
 template <int dim>
 std::string
 FE_Q<dim>::get_name () const
@@ -221,8 +274,62 @@ FE_Q<dim>::get_name () const
                                   // have to be kept in synch
   
   std::ostringstream namebuf;  
-  namebuf << "FE_Q<" << dim << ">(" << this->degree << ")";
+  bool type = true;  
+  const unsigned int n_points = this->degree +1;
+  std::vector<double> points(n_points);
+  const unsigned int dofs_per_cell = this->dofs_per_cell;
+  const std::vector<Point<dim> > &unit_support_points = this->unit_support_points;
+  unsigned int index = 0;
+
+                                  // Decode the support points
+                                  // in one coordinate direction.
+  for (unsigned int j=0;j<dofs_per_cell;j++)
+    {
+      if ((dim>1) ? (unit_support_points[j](1)==0 && 
+          ((dim>2) ? unit_support_points[j](2)==0: true)) : true)
+       {
+         if (index == 0)
+           points[index] = unit_support_points[j](0);
+          else if (index == 1)
+           points[n_points-1] = unit_support_points[j](0);
+          else
+           points[index-1] = unit_support_points[j](0);
 
+         index++;
+       }
+    }
+  Assert (index == n_points,
+         ExcMessage ("Could not decode support points in one coordinate direction."));
+
+                                 // Check whether the support
+                                 // points are equidistant.
+  for(unsigned int j=0;j<n_points;j++)
+    if (std::fabs(points[j] - (double)j/this->degree) > 1e-15)
+      {
+       type = false;
+       break;
+      }
+
+  if (type == true)    
+    namebuf << "FE_Q<" << dim << ">(" << this->degree << ")";
+  else
+    {
+
+                                 // Check whether the support
+                                 // points come from QGaussLobatto.
+      const QGaussLobatto<1> points_gl(n_points);
+      type = true;
+      for(unsigned int j=0;j<n_points;j++)
+       if (points[j] != points_gl.point(j)(0))
+         {
+           type = false;
+           break;
+         }
+      if(type == true)
+       namebuf << "FE_Q<" << dim << ">(QGaussLobatto(" << this->degree+1 << "))";
+      else
+       namebuf << "FE_Q<" << dim << ">(QUnknownNodes(" << this->degree << "))";
+    }
   return namebuf.str();
 }
 
@@ -241,7 +348,7 @@ template <int dim>
 void
 FE_Q<dim>::
 get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
-                         FullMatrix<double>           &interpolation_matrix) const
+                         FullMatrix<double>       &interpolation_matrix) const
 {
                                   // this is only implemented, if the
                                   // source FE is also a
@@ -264,9 +371,6 @@ get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
   const FE_Q<dim> &source_fe
     = dynamic_cast<const FE_Q<dim>&>(x_source_fe);
 
-  const std::vector<unsigned int> &index_map=
-    this->poly_space.get_numbering();
-  
                                   // compute the interpolation
                                   // matrices in much the same way as
                                   // we do for the embedding matrices
@@ -279,12 +383,10 @@ get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
                          source_fe.dofs_per_cell);
   for (unsigned int j=0; j<this->dofs_per_cell; ++j)
     {
-                                   // generate a point on this
+                                   // read in a point on this
                                    // cell and evaluate the
                                    // shape functions there
-      const Point<dim>
-       p = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell,
-                                             dealii::internal::int2type<dim>());
+      const Point<dim> p = this->unit_support_points[j];
       for (unsigned int i=0; i<this->dofs_per_cell; ++i)
         cell_interpolation(j,i) = this->poly_space.compute_value (i, p);
 
@@ -293,7 +395,7 @@ get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
     }
 
                                    // then compute the
-                                   // interpolation matrix matrix
+                                   // interpolation matrix
                                    // for this coordinate
                                    // direction
   cell_interpolation.gauss_jordan ();
@@ -337,6 +439,7 @@ get_face_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/,
 }
 
 
+
 template <>
 void
 FE_Q<1>::
@@ -424,15 +527,15 @@ get_face_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
       const Point<dim> &p = face_quadrature.point (i);
 
       for (unsigned int j=0; j<this->dofs_per_face; ++j)
-       { 
+       {
          double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
 
-                                          // Correct the interpolated
-                                          // value. I.e. if it is close
-                                          // to 1 or 0, make it exactly
-                                          // 1 or 0. Unfortunately, this
-                                          // is required to avoid problems
-                                          // with higher order elements.
+                                  // Correct the interpolated
+                                  // value. I.e. if it is close
+                                  // to 1 or 0, make it exactly
+                                  // 1 or 0. Unfortunately, this
+                                  // is required to avoid problems
+                                  // with higher order elements.
          if (fabs (matrix_entry - 1.0) < eps)
            matrix_entry = 1.0;
          if (fabs (matrix_entry) < eps)
@@ -460,12 +563,13 @@ get_face_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
 }
 
 
+
 template <int dim>
 void
 FE_Q<dim>::
 get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
-                                 const unsigned int subface,
-                                 FullMatrix<double>           &interpolation_matrix) const
+                                 const unsigned int        subface,
+                                 FullMatrix<double>       &interpolation_matrix) const
 {
                                   // this is only implemented, if the
                                   // source FE is also a
@@ -747,6 +851,31 @@ void FE_Q<dim>::initialize_unit_support_points ()
 }
 
 
+
+template <int dim>
+void FE_Q<dim>::initialize_unit_support_points (const Quadrature<1> &points)
+{
+                                  // number of points: (degree+1)^dim
+  unsigned int n = this->degree+1;
+  for (unsigned int i=1; i<dim; ++i)
+    n *= this->degree+1;
+  
+  this->unit_support_points.resize(n);
+
+  const std::vector<unsigned int> &index_map_inverse=
+    this->poly_space.get_numbering_inverse();
+       
+  Quadrature<dim> support_quadrature(points);
+
+  Point<dim> p;
+  
+  for (unsigned int k=0;k<n ;k++)
+    {
+      this->unit_support_points[index_map_inverse[k]] = support_quadrature.point(k);
+    }
+}
+
+
 #if deal_II_dimension == 1
 
 template <>
@@ -755,6 +884,12 @@ void FE_Q<1>::initialize_unit_face_support_points ()
                                   // no faces in 1d, so nothing to do
 }
 
+template <>
+void FE_Q<1>::initialize_unit_face_support_points (const Quadrature<1> &/*points*/)
+{
+                                  // no faces in 1d, so nothing to do
+}
+
 #endif
 
 
@@ -793,6 +928,42 @@ void FE_Q<dim>::initialize_unit_face_support_points ()
 
 
 
+template <int dim>
+void FE_Q<dim>::initialize_unit_face_support_points (const Quadrature<1> &points)
+{
+  const unsigned int codim = dim-1;
+  
+                                  // number of points: (degree+1)^codim
+  unsigned int n = this->degree+1;
+  for (unsigned int i=1; i<codim; ++i)
+    n *= this->degree+1;
+  
+  this->unit_face_support_points.resize(n);
+
+  const std::vector< Point<1> > edge = points.get_points();
+
+  const std::vector<unsigned int> &face_index_map_inverse=
+    FE_Q_Helper::invert_numbering(face_index_map);
+  
+  Point<codim> p;
+  
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=this->degree; ++ix)
+       {
+         p(0) = edge[ix](0);
+         if (codim>1)
+           p(1) = edge[iy](0);
+         if (codim>2)
+           p(2) = edge[iz](0);
+         
+         this->unit_face_support_points[face_index_map_inverse[k++]] = p;
+       }
+}
+
+
+
 template <int dim>
 void
 FE_Q<dim>::initialize_quad_dof_index_permutation ()
@@ -801,7 +972,6 @@ FE_Q<dim>::initialize_quad_dof_index_permutation ()
 }
 
 
-
 #if deal_II_dimension == 3
 
 template <>
@@ -885,6 +1055,7 @@ FE_Q<dim>::get_dpo_vector(const unsigned int deg)
 }
 
 
+
 template <int dim>
 std::vector<unsigned int>
 FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
@@ -907,22 +1078,35 @@ FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
 
 #endif
 
+
+
+template <int dim>
+void
+FE_Q<dim>::initialize_constraints ()
+{
+  QTrapez<1> trapez;
+  QIterated<1> points (trapez,this->degree);
+  initialize_constraints (points);
+}
+
+
 #if deal_II_dimension == 1
 
 template <>
 void
-FE_Q<1>::initialize_constraints ()
+FE_Q<1>::initialize_constraints (const Quadrature<1> &/*points*/)
 {
                                   // no constraints in 1d
 }
 
 #endif
 
+
 #if deal_II_dimension == 2
 
 template <>
 void
-FE_Q<2>::initialize_constraints ()
+FE_Q<2>::initialize_constraints (const Quadrature<1> &points)
 {
   const unsigned int dim = 2;
 
@@ -1022,7 +1206,7 @@ FE_Q<2>::initialize_constraints ()
                                    // destination (child) and source (mother)
                                    // dofs.
   const std::vector<Polynomials::Polynomial<double> > polynomials=
-    Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
+    Polynomials::Lagrange::generate_complete_basis(points.get_points());
 
   this->interface_constraints
     .TableBase<2,double>::reinit (this->interface_constraints_size());
@@ -1033,16 +1217,16 @@ FE_Q<2>::initialize_constraints ()
        this->interface_constraints(i,j) = 
          polynomials[face_index_map[j]].value (constraint_points[i](0));
                    
-                                           // if the value is small up
-                                           // to round-off, then
-                                           // simply set it to zero to
-                                           // avoid unwanted fill-in
-                                           // of the constraint
-                                           // matrices (which would
-                                           // then increase the number
-                                           // of other DoFs a
-                                           // constrained DoF would
-                                           // couple to)
+                                   // if the value is small up
+                                   // to round-off, then
+                                  // simply set it to zero to
+                                   // avoid unwanted fill-in
+                                   // of the constraint
+                                   // matrices (which would
+                                   // then increase the number
+                                   // of other DoFs a
+                                   // constrained DoF would
+                                   // couple to)
           if (std::fabs(this->interface_constraints(i,j)) < 1e-14)
             this->interface_constraints(i,j) = 0;
       }
@@ -1051,9 +1235,10 @@ FE_Q<2>::initialize_constraints ()
 #endif
 
 #if deal_II_dimension == 3
+
 template <>
 void
-FE_Q<3>::initialize_constraints ()
+FE_Q<3>::initialize_constraints (const Quadrature<1> &points)
 {
   const unsigned int dim = 3;
 
@@ -1148,9 +1333,11 @@ FE_Q<3>::initialize_constraints ()
                                    // Now construct relation between
                                    // destination (child) and source (mother)
                                    // dofs.
-  const unsigned int pnts=(this->degree+1)*(this->degree+1);  
+  const unsigned int pnts=(this->degree+1)*(this->degree+1); 
   const std::vector<Polynomials::Polynomial<double> > polynomial_basis=
-    Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
+    Polynomials::Lagrange::generate_complete_basis(points.get_points()); 
+    //const std::vector<Polynomials::Polynomial<double> > polynomial_basis=
+    //Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
  
   const TensorProductPolynomials<dim-1> face_polynomials(polynomial_basis);
 
index 83f3e21908995923dac176d262651f2c4bfe0af5..42dd5e187204f2bb0d21ade4320e90955545ffbc 100644 (file)
@@ -92,6 +92,8 @@ namespace
       = FEFactoryPointer(new FETools::FEFactory<FE_DGPMonomial<deal_II_dimension> >);
     default_map["FE_DGQ"]
       = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<deal_II_dimension> >);
+    default_map["FE_DGQArbitraryNodes"]
+      = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<deal_II_dimension> >);
     default_map["FE_Q"]
       = FEFactoryPointer(new FETools::FEFactory<FE_Q<deal_II_dimension> >);
 
@@ -1751,14 +1753,37 @@ FETools::get_fe_from_name_aux (std::string &name)
        AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(),
                      ExcInvalidFEName(name));
                                         // Now, just the (degree)
+                                        // or (Quadrature<1>(degree+1))
                                         // part should be left.
        if (name.size() == 0 || name[0] != '(')
          throw (std::string("Invalid first character in ") + name);
        name.erase(0,1);
-       const std::pair<int,unsigned int> tmp
-         = Utilities::get_integer_at_position (name, 0);
-       name.erase(0, tmp.second+1);
-       return fe_name_map.find(name_part)->second->get(tmp.first);
+       if (name[0] != 'Q')
+         {
+           const std::pair<int,unsigned int> tmp
+             = Utilities::get_integer_at_position (name, 0);
+           name.erase(0, tmp.second+1);
+           return fe_name_map.find(name_part)->second->get(tmp.first);
+         }
+        else
+          {
+            unsigned int position = name.find('(');
+            const std::string quadrature_name(name, 0, position-1);
+            name.erase(0,position);
+            if (quadrature_name.compare("QGaussLobatto") == 0)         
+              {
+                const std::pair<int,unsigned int> tmp
+                  = Utilities::get_integer_at_position (name, 0);
+                name.erase(0, tmp.second+1);
+//TODO: Implement a get function taking Quadrature<1> in fe_tools.h.
+                //return fe_name_map.find(name_part)->second->get(QGaussLobatto<1>(tmp.first));
+                AssertThrow (false, ExcNotImplemented());
+              }
+            else 
+              {
+                AssertThrow (false,ExcNotImplemented());
+              }
+          }
       }
   
     
index 72d326ba7c435e120967cf07013be9749e5878d9..6d2d659e29975fecb623d8bef8cfce07dd1d86b3 100644 (file)
@@ -95,6 +95,9 @@ alphabetical order):
 <li><em>Benjamin Shelton Kirk:</em>
     Tecplot output.
 
+<li><em>Katharina Kormann:</em>
+    Support for arbitrary nodes in FE_Q.
+
 <li><em>Martin Kronbichler:</em>
     Parts of step-22 and step-31 tutorial programs, interfaces to Trilinos (partly), and a few enhancements.
 
index 18222d823153828bd7f4d0a88bb50cb3e20cb1fe..393600768981c7ecab0efc920492403f4911a24d 100644 (file)
@@ -278,6 +278,14 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>
+  New: The class FE_Q can now alternatively be constructed based on 
+  support points from a given one-dimensional quadrature rule.
+  <br>
+  (Katharina Kormann, Martin Kronbichler, 2008/09/07)
+  </p> 
+
   <li>
   <p>
   Fixed: Using the ConstraintMatrix class, when a degree of freedom was

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.