]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
do some cleaning in the RT elements
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jan 2009 20:40:08 +0000 (20:40 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jan 2009 20:40:08 +0000 (20:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@18164 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_raviart_thomas.h
deal.II/deal.II/source/fe/fe_poly_tensor.cc
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc

index fdb283259de50de4a545d969af40e366e0ae466e..fde8131c31288db2c87d2b85b1c9a7c276ae2f89 100644 (file)
@@ -112,34 +112,6 @@ class FE_RaviartThomas
     virtual std::string get_name () const;
 
 
-                                    /**
-                                     * Number of base elements in a
-                                     * mixed discretization. Here,
-                                     * this is of course equal to
-                                     * one.
-                                     */
-    virtual unsigned int n_base_elements () const;
-    
-                                    /**
-                                     * Access to base element
-                                     * objects. Since this element is
-                                     * atomic, <tt>base_element(0)</tt> is
-                                     * @p this, and all other
-                                     * indices throw an error.
-                                     */
-    virtual const FiniteElement<dim> &
-    base_element (const unsigned int index) const;
-
-                                     /**
-                                      * Multiplicity of base element
-                                      * @p index. Since this is an
-                                      * atomic element,
-                                      * <tt>element_multiplicity(0)</tt>
-                                      * returns one, and all other
-                                      * indices will throw an error.
-                                      */
-    virtual unsigned int element_multiplicity (const unsigned int index) const;
-    
                                     /**
                                      * Check whether a shape function
                                      * may be non-zero on a face.
@@ -164,17 +136,6 @@ class FE_RaviartThomas
     virtual FiniteElement<dim> * clone() const;
     
   private:
-                                    /**
-                                     * The order of the
-                                     * Raviart-Thomas element. The
-                                     * lowest order elements are
-                                     * usually referred to as RT0,
-                                     * even though their shape
-                                     * functions are piecewise
-                                     * linears.
-                                     */  
-    const unsigned int rt_order;
-
                                     /**
                                      * Only for internal use. Its
                                      * full name is
@@ -213,45 +174,6 @@ class FE_RaviartThomas
                                      */
     void initialize_restriction ();
     
-                                    /**
-                                     * Given a set of flags indicating
-                                     * what quantities are requested
-                                     * from a @p FEValues object,
-                                     * return which of these can be
-                                     * precomputed once and for
-                                     * all. Often, the values of
-                                     * shape function at quadrature
-                                     * points can be precomputed, for
-                                     * example, in which case the
-                                     * return value of this function
-                                     * would be the logical and of
-                                     * the input @p flags and
-                                     * @p update_values.
-                                     *
-                                     * For the present kind of finite
-                                     * element, this is exactly the
-                                     * case.
-                                     */
-    virtual UpdateFlags update_once (const UpdateFlags flags) const;
-  
-                                    /**
-                                     * This is the opposite to the
-                                     * above function: given a set of
-                                     * flags indicating what we want
-                                     * to know, return which of these
-                                     * need to be computed each time
-                                     * we visit a new cell.
-                                     *
-                                     * If for the computation of one
-                                     * quantity something else is
-                                     * also required (for example, we
-                                     * often need the covariant
-                                     * transformation when gradients
-                                     * need to be computed), include
-                                     * this in the result as well.
-                                     */
-    virtual UpdateFlags update_each (const UpdateFlags flags) const;
-    
                                     /**
                                      * Fields of cell-independent data.
                                      *
@@ -442,8 +364,6 @@ class FE_RaviartThomasNodal
     virtual FiniteElementDomination::Domination
     compare_for_face_domination (const FiniteElement<dim> &fe_other) const;
 
-    virtual UpdateFlags update_once (const UpdateFlags flags) const;
-    virtual UpdateFlags update_each (const UpdateFlags flags) const;
   private:
                                     /**
                                      * Only for internal use. Its
index 1ce943a7c752265da2ac15e22a260677eed8aae4..9e1dc328fc2e68dfc800b71af97d5994f5066ec2 100644 (file)
@@ -860,19 +860,30 @@ template <class POLY, int dim, int spacedim>
 UpdateFlags
 FE_PolyTensor<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const
 {
-  const bool values_once = (mapping_type == mapping_none);
-  
   UpdateFlags out = update_default;
 
-  if (!values_once && (flags & update_values))
-    out |= update_values             | update_covariant_transformation;
-  if (flags & update_gradients)
-    out |= update_gradients          | update_covariant_transformation;
-  if (flags & update_hessians)
-    out |= update_hessians | update_covariant_transformation;
-  if (flags & update_cell_normal_vectors)
-    out |= update_cell_normal_vectors | update_JxW_values;
-
+  switch (mapping_type)
+    {
+      case mapping_piola:
+      {
+       if (flags & update_values)
+         out |= update_values | update_piola;
+       
+       if (flags & update_gradients)
+         out |= update_gradients | update_piola | update_covariant_transformation;
+       
+       if (flags & update_hessians)
+         out |= update_hessians | update_piola | update_covariant_transformation;
+       
+       break;
+      }
+      
+      default:
+      {
+       Assert (false, ExcNotImplemented());
+      }
+    }
+  
   return out;
 }
 
index cce73d0fcc699558abeb041d23803060ba07da87..cf158ff4b2460b178ea57ed75ea9c445e7470c20 100644 (file)
@@ -35,6 +35,7 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+
 template <int dim>
 FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
                :
@@ -44,8 +45,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
                                         dim, deg+1, FiniteElementData<dim>::Hdiv, 1),
                  std::vector<bool>(PolynomialsRaviartThomas<dim>::compute_n_pols(deg), true),
                  std::vector<std::vector<bool> >(PolynomialsRaviartThomas<dim>::compute_n_pols(deg),
-                                                 std::vector<bool>(dim,true))),
-               rt_order(deg)
+                                                 std::vector<bool>(dim,true)))
 {
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
@@ -60,6 +60,12 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
                                   //matrix, generating the correct
                                   //basis functions from the raw
                                   //ones.
+  
+                                  // We use an auxiliary matrix in
+                                  // this function. Therefore,
+                                  // inverse_node_matrix is still
+                                  // empty and shape_value_component
+                                  // returns the 'raw' shape values.
   FullMatrix<double> M(n_dofs, n_dofs);
   FETools::compute_node_matrix(M, *this);
   this->inverse_node_matrix.reinit(n_dofs, n_dofs);
@@ -82,9 +88,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i=0; i<GeometryInfo<dim>::max_children_per_face; ++i)
     face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
-  FETools::compute_face_embedding_matrices<dim,double>(*this, 
-                                                      face_embeddings, 
-                                                      0, 0);
+  FETools::compute_face_embedding_matrices<dim,double>(*this, face_embeddings, 0, 0);
   this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face,
                                     this->dofs_per_face);
   unsigned int target_row=0;
@@ -110,14 +114,17 @@ FE_RaviartThomas<dim>::get_name () const
                                   // this function returns, so they
                                   // have to be kept in synch
 
+                                  // note that this->degree is the maximal
+                                  // polynomial degree and is thus one higher
+                                  // than the argument given to the
+                                  // constructor
   std::ostringstream namebuf;  
-  namebuf << "FE_RaviartThomas<" << dim << ">(" << rt_order << ")";
+  namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree-1 << ")";
 
   return namebuf.str();
 }
 
 
-
 template <int dim>
 FiniteElement<dim> *
 FE_RaviartThomas<dim>::clone() const
@@ -288,7 +295,7 @@ FE_RaviartThomas<dim>::initialize_restriction()
 {
   const unsigned int iso=RefinementCase<dim>::isotropic_refinement-1;
 
-  QGauss<dim-1> q_base (rt_order+1);
+  QGauss<dim-1> q_base (this->degree);
   const unsigned int n_face_points = q_base.size();
                                   // First, compute interpolation on
                                   // subfaces
@@ -354,7 +361,7 @@ FE_RaviartThomas<dim>::initialize_restriction()
        }
     }
   
-  if (rt_order==0) return;
+  if (this->degree == 1) return;
   
                                   // Create Legendre basis for the
                                   // space D_xi Q_k. Here, we cannot
@@ -364,13 +371,13 @@ FE_RaviartThomas<dim>::initialize_restriction()
     {
       std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
       for (unsigned int d=0;d<dim;++d)
-       poly[d] = Polynomials::Legendre::generate_complete_basis(rt_order);
-      poly[dd] = Polynomials::Legendre::generate_complete_basis(rt_order-1);
+       poly[d] = Polynomials::Legendre::generate_complete_basis(this->degree-1);
+      poly[dd] = Polynomials::Legendre::generate_complete_basis(this->degree-2);
 
       polynomials[dd] = new AnisotropicPolynomials<dim>(poly);
     }
   
-  QGauss<dim> q_cell(rt_order+1);
+  QGauss<dim> q_cell(this->degree);
   const unsigned int start_cell_dofs
     = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
 
@@ -444,72 +451,11 @@ FE_RaviartThomas<dim>::get_dpo_vector (const unsigned int rt_order)
 
 
 
-template <int dim>
-UpdateFlags
-FE_RaviartThomas<dim>::update_once (const UpdateFlags) const
-{
-                                  // even the values have to be
-                                  // computed on the real cell, so
-                                  // nothing can be done in advance
-  return update_default;
-}
-
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
-{
-  UpdateFlags out = update_default;
-
-  if (flags & update_values)
-    out |= update_values | update_piola;
-  
-  if (flags & update_gradients)
-    out |= update_gradients | update_piola | update_covariant_transformation;
-  
-  if (flags & update_hessians)
-    out |= update_hessians | update_piola | update_covariant_transformation;
-  
-  return out;
-}
-
 //---------------------------------------------------------------------------
 // Data field initialization
 //---------------------------------------------------------------------------
 
 
-
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::n_base_elements () const
-{
-  return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_RaviartThomas<dim>::base_element (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::element_multiplicity (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return 1;
-}
-
-
-
 template <int dim>
 bool
 FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
@@ -523,9 +469,9 @@ FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
                                   // Return computed values if we
                                   // know them easily. Otherwise, it
                                   // is always safe to return true.
-  switch (rt_order)
+  switch (this->degree)
     {
-      case 0:
+      case 1:
       {
         switch (dim)
           {
index 9c0087f7cd4d7bc1b7129a7b3cbb7c8b8b0b0681..8eb0fa196c73eb3504fc8144bf5133a69271d74d 100644 (file)
@@ -25,6 +25,7 @@
 
 #include <sstream>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -43,14 +44,12 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
   const unsigned int n_dofs = this->dofs_per_cell;
   
   this->mapping_type = mapping_raviart_thomas;
-                                  // These must be done first, since
-                                  // they change the evaluation of
-                                  // basis functions
-
-                                  // Set up the generalized support
-                                  // points
-  initialize_support_points (deg);
-                                  //Now compute the inverse node
+                                  // First, initialize the
+                                  // generalized support points and
+                                  // quadrature weights, since they
+                                  // are required for interpolation.
+  initialize_support_points(deg);
+                                  // Now compute the inverse node
                                   //matrix, generating the correct
                                   //basis functions from the raw
                                   //ones.
@@ -131,6 +130,147 @@ FE_RaviartThomasNodal<dim>::clone() const
 }
 
 
+//---------------------------------------------------------------------------
+// Auxiliary and internal functions
+//---------------------------------------------------------------------------
+
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::initialize_support_points (const unsigned int deg)
+{
+  this->generalized_support_points.resize (this->dofs_per_cell);
+  this->generalized_face_support_points.resize (this->dofs_per_face);
+
+                                  // Number of the point being entered
+  unsigned int current = 0;
+
+                                  // On the faces, we choose as many
+                                  // Gauss points as necessary to
+                                  // determine the normal component
+                                  // uniquely. This is the deg of
+                                  // the Raviart-Thomas element plus
+                                  // one.
+  if (dim>1)
+    {
+      QGauss<dim-1> face_points (deg+1);
+      Assert (face_points.size() == this->dofs_per_face,
+             ExcInternalError());
+      for (unsigned int k=0;k<this->dofs_per_face;++k)
+       this->generalized_face_support_points[k] = face_points.point(k);
+      Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
+      for (unsigned int k=0;
+          k<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
+          ++k)
+       this->generalized_support_points[k] = faces.point(k+QProjector<dim>
+                                                         ::DataSetDescriptor::face(0,
+                                                                                   true,
+                                                                                   false,
+                                                                                   false,
+                                                                                   this->dofs_per_face));
+
+      current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
+    }
+  
+  if (deg==0) return;
+                                  // In the interior, we need
+                                  // anisotropic Gauss quadratures,
+                                  // different for each direction.
+  QGauss<1> high(deg+1);
+  QGauss<1> low(deg);
+
+  for (unsigned int d=0;d<dim;++d)
+    {
+      QAnisotropic<dim>* quadrature;
+      if (dim == 1) quadrature = new QAnisotropic<dim>(high);
+      if (dim == 2) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
+                                                      ((d==1) ? low : high));
+      if (dim == 3) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
+                                                      ((d==1) ? low : high),
+                                                      ((d==2) ? low : high));
+      Assert(dim<=3, ExcNotImplemented());
+      
+      for (unsigned int k=0;k<quadrature->size();++k)
+       this->generalized_support_points[current++] = quadrature->point(k);
+      delete quadrature;
+    }
+  Assert (current == this->dofs_per_cell, ExcInternalError());
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<unsigned int>
+FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg)
+{
+  std::vector<unsigned int> dpo(2);
+  dpo[0] = 1;
+  dpo[1] = deg;
+  return dpo;
+}
+
+#endif
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_RaviartThomasNodal<dim>::get_dpo_vector (const unsigned int deg)
+{
+                                   // the element is face-based and we have
+                                   // (deg+1)^(dim-1) DoFs per face
+  unsigned int dofs_per_face = 1;
+  for (unsigned int d=1; d<dim; ++d)
+    dofs_per_face *= deg+1;
+
+                                   // and then there are interior dofs
+  const unsigned int
+    interior_dofs = dim*deg*dofs_per_face;
+  
+  std::vector<unsigned int> dpo(dim+1);
+  dpo[dim-1] = dofs_per_face;
+  dpo[dim]   = interior_dofs;
+  
+  return dpo;
+}
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<bool>
+FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int)
+{
+  Assert (false, ExcImpossibleInDim(1));
+  return std::vector<bool>();
+}
+
+#endif
+
+
+template <int dim>
+std::vector<bool>
+FE_RaviartThomasNodal<dim>::get_ria_vector (const unsigned int deg)
+{
+  const unsigned int dofs_per_cell = PolynomialsRaviartThomas<dim>::compute_n_pols(deg);
+  unsigned int dofs_per_face = deg+1;
+  for (unsigned int d=2;d<dim;++d)
+    dofs_per_face *= deg+1;
+                                  // all face dofs need to be
+                                  // non-additive, since they have
+                                  // continuity requirements.
+                                  // however, the interior dofs are
+                                  // made additive
+  std::vector<bool> ret_val(dofs_per_cell,false);
+  for (unsigned int i=GeometryInfo<dim>::faces_per_cell*dofs_per_face;
+       i < dofs_per_cell; ++i)
+    ret_val[i] = true;
+
+  return ret_val;
+}
+
 
 template <int dim>
 void
@@ -634,169 +774,6 @@ get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
 
 
 
-#if deal_II_dimension == 1
-
-template <>
-std::vector<unsigned int>
-FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg)
-{
-  std::vector<unsigned int> dpo(2);
-  dpo[0] = 1;
-  dpo[1] = deg;
-  return dpo;
-}
-
-#endif
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_RaviartThomasNodal<dim>::get_dpo_vector (const unsigned int deg)
-{
-                                   // the element is face-based and we have
-                                   // (deg+1)^(dim-1) DoFs per face
-  unsigned int dofs_per_face = 1;
-  for (unsigned int d=1; d<dim; ++d)
-    dofs_per_face *= deg+1;
-
-                                   // and then there are interior dofs
-  const unsigned int
-    interior_dofs = dim*deg*dofs_per_face;
-  
-  std::vector<unsigned int> dpo(dim+1);
-  dpo[dim-1] = dofs_per_face;
-  dpo[dim]   = interior_dofs;
-  
-  return dpo;
-}
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-std::vector<bool>
-FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int)
-{
-  Assert (false, ExcImpossibleInDim(1));
-  return std::vector<bool>();
-}
-
-#endif
-
-
-template <int dim>
-std::vector<bool>
-FE_RaviartThomasNodal<dim>::get_ria_vector (const unsigned int deg)
-{
-  const unsigned int dofs_per_cell = PolynomialsRaviartThomas<dim>::compute_n_pols(deg);
-  unsigned int dofs_per_face = deg+1;
-  for (unsigned int d=2;d<dim;++d)
-    dofs_per_face *= deg+1;
-                                  // all face dofs need to be
-                                  // non-additive, since they have
-                                  // continuity requirements.
-                                  // however, the interior dofs are
-                                  // made additive
-  std::vector<bool> ret_val(dofs_per_cell,false);
-  for (unsigned int i=GeometryInfo<dim>::faces_per_cell*dofs_per_face;
-       i < dofs_per_cell; ++i)
-    ret_val[i] = true;
-
-  return ret_val;
-}
-
-
-template <int dim>
-void
-FE_RaviartThomasNodal<dim>::initialize_support_points (const unsigned int deg)
-{
-  this->generalized_support_points.resize (this->dofs_per_cell);
-  this->generalized_face_support_points.resize (this->dofs_per_face);
-
-                                  // Number of the point being entered
-  unsigned int current = 0;
-
-                                  // On the faces, we choose as many
-                                  // Gauss points as necessary to
-                                  // determine the normal component
-                                  // uniquely. This is the deg of
-                                  // the Raviart-Thomas element plus
-                                  // one.
-  if (dim>1)
-    {
-      QGauss<dim-1> face_points (deg+1);
-      Assert (face_points.size() == this->dofs_per_face,
-             ExcInternalError());
-      for (unsigned int k=0;k<this->dofs_per_face;++k)
-       this->generalized_face_support_points[k] = face_points.point(k);
-      Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
-      for (unsigned int k=0;
-          k<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
-          ++k)
-       this->generalized_support_points[k] = faces.point(k+QProjector<dim>
-                                                         ::DataSetDescriptor::face(0,
-                                                                                   true,
-                                                                                   false,
-                                                                                   false,
-                                                                                   this->dofs_per_face));
-
-      current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
-    }
-  
-  if (deg==0) return;
-                                  // In the interior, we need
-                                  // anisotropic Gauss quadratures,
-                                  // different for each direction.
-  QGauss<1> high(deg+1);
-  QGauss<1> low(deg);
-
-  for (unsigned int d=0;d<dim;++d)
-    {
-      QAnisotropic<dim>* quadrature;
-      if (dim == 1) quadrature = new QAnisotropic<dim>(high);
-      if (dim == 2) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
-                                                      ((d==1) ? low : high));
-      if (dim == 3) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
-                                                      ((d==1) ? low : high),
-                                                      ((d==2) ? low : high));
-      Assert(dim<=3, ExcNotImplemented());
-      
-      for (unsigned int k=0;k<quadrature->size();++k)
-       this->generalized_support_points[current++] = quadrature->point(k);
-      delete quadrature;
-    }
-  Assert (current == this->dofs_per_cell, ExcInternalError());
-}
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomasNodal<dim>::update_once (const UpdateFlags) const
-{
-  return update_default;
-}
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomasNodal<dim>::update_each (const UpdateFlags flags) const
-{
-  UpdateFlags out = update_default;
-  
-  if (flags & update_values)
-    out |= update_values | update_piola;
-  
-  if (flags & update_gradients)
-    out |= update_gradients | update_piola | update_covariant_transformation;
-  
-  if (flags & update_hessians)
-    out |= update_hessians | update_piola | update_covariant_transformation;
-  
-  return out;
-}
-
-
 template class FE_RaviartThomasNodal<deal_II_dimension>;
 
 DEAL_II_NAMESPACE_CLOSE

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.