]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new RT compiles and links
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 24 May 2005 10:20:27 +0000 (10:20 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 24 May 2005 10:20:27 +0000 (10:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@10719 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_poly_tensor.h
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_nodal.cc [new file with mode: 0644]

index f18a0061f55565631558b82f632e918cc15f190e..4c8559023bc42083fc8ef6419f3edfc0a14cb4cf 100644 (file)
@@ -63,6 +63,11 @@ class FE_PolyTensor : public FiniteElement<dim>
   public:
                                     /**
                                      * Constructor.
+                                     *
+                                     * @arg @c degree: constructor
+                                     * argument for poly. May be
+                                     * different from @p
+                                     * fe_data.degree.
                                      */
     FE_PolyTensor (unsigned int degree,
                   const FiniteElementData<dim> &fe_data,
@@ -134,6 +139,28 @@ class FE_PolyTensor : public FiniteElement<dim>
                                       */
     virtual unsigned int element_multiplicity (const unsigned int index) const;
 
+                                    /**
+                                     * Given <tt>flags</tt>,
+                                     * determines the values which
+                                     * must be computed only for the
+                                     * reference cell. Make sure,
+                                     * that #mapping_type is set by
+                                     * the derived class, such that
+                                     * this function can operate
+                                     * correctly.
+                                     */
+    virtual UpdateFlags update_once (const UpdateFlags flags) const;
+                                    /**
+                                     * Given <tt>flags</tt>,
+                                     * determines the values which
+                                     * must be computed in each cell
+                                     * cell. Make sure, that
+                                     * #mapping_type is set by the
+                                     * derived class, such that this
+                                     * function can operate
+                                     * correctly.
+                                     */
+    virtual UpdateFlags update_each (const UpdateFlags flags) const;
     
   protected:
                                     /**
@@ -150,11 +177,36 @@ class FE_PolyTensor : public FiniteElement<dim>
                                      * possible to avoid the mapping.
                                      */
     enum MappingType {
-                                          /// Shape functions do not depend on actual mesh cell
+                                          /**
+                                           * No mapping has been
+                                           * selected, throw an error
+                                           * if needed.
+                                           */
+         no_mapping,
+                                          /**
+                                           * Shape functions do not
+                                           * depend on actual mesh
+                                           * cell
+                                           */
          independent,
-                                          /// Shape functions are transformed covariant.
+                                          /**
+                                           * Shape functions do not
+                                           * depend on actual mesh
+                                           * cell. The mapping class
+                                           * must be
+                                           * MappingCartesian.
+                                           */
+         independent_on_cartesian,
+                                          /**
+                                           * Shape functions are
+                                           * transformed covariant.
+                                           */ 
          covariant,
-                                          /// Shape functions are transformed contravariant.
+                                          /**
+                                           * Shape functions are
+                                           * transformed
+                                           * contravariant.
+                                           */
          contravariant
     };
 
@@ -240,17 +292,13 @@ class FE_PolyTensor : public FiniteElement<dim>
        std::vector<std::vector<Tensor<2,dim> > > shape_grads;
     };
     
-                                    /**
-                                     * Degree of the polynomials.
-                                     */  
-    unsigned int degree;
-
                                      /**
                                       * The polynomial space. Its type
                                       * is given by the template
                                       * parameter POLY.
                                       */    
     POLY poly_space;
+    
                                     /**
                                      * The inverse of the matrix
                                      * <i>a<sub>ij</sub></i> of node
index d9b52ead798037b64a2caabd65d2da97ba662d90..3a8cf09a528c3262fa74eb71d2fd49d4bf000451 100644 (file)
 #define __deal2__fe_raviart_thomas_h
 
 #include <base/config.h>
+#include <base/polynomials_raviart_thomas.h>
 #include <base/tensor_product_polynomials.h>
 #include <grid/geometry_info.h>
 #include <fe/fe.h>
+#include <fe/fe_poly_tensor.h>
 
 template <int dim> class MappingQ;
 
@@ -258,16 +260,12 @@ class FE_RaviartThomas : public FiniteElement<dim>
     
                                     /**
                                      * Check whether a shape function
-                                     * is non-zero on a face.
+                                     * may be non-zero on a face.
                                      *
                                      * Right now, this is only
                                      * implemented for RT0 in
                                      * 1D. Otherwise, returns always
                                      * @p true.
-                                     *
-                                     * Implementation of the
-                                     * interface in
-                                     * FiniteElement
                                      */
     virtual bool has_support_on_face (const unsigned int shape_index,
                                      const unsigned int face_index) const;
@@ -291,31 +289,15 @@ class FE_RaviartThomas : public FiniteElement<dim>
     DeclException0 (ExcNotUsefulInThisDimension);
     
   protected:    
-                                    /**
-                                     * @p clone function instead of
-                                     * a copy constructor.
-                                     *
-                                     * This function is needed by the
-                                     * constructors of @p FESystem.
-                                     */
+    
     virtual FiniteElement<dim> * clone() const;
   
-                                    /**
-                                     * Prepare internal data
-                                     * structures and fill in values
-                                     * independent of the cell.
-                                     */
     virtual
     typename Mapping<dim>::InternalDataBase *
     get_data (const UpdateFlags,
              const Mapping<dim>& mapping,
              const Quadrature<dim>& quadrature) const ;
 
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * FiniteElement.
-                                     */
     virtual void
     fill_fe_values (const Mapping<dim> &mapping,
                    const typename Triangulation<dim>::cell_iterator &cell,
@@ -324,11 +306,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
                    typename Mapping<dim>::InternalDataBase      &fe_internal,
                    FEValuesData<dim>& data) const;
     
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * FiniteElement.
-                                     */
     virtual void
     fill_fe_face_values (const Mapping<dim> &mapping,
                         const typename Triangulation<dim>::cell_iterator &cell,
@@ -338,11 +315,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
                         typename Mapping<dim>::InternalDataBase      &fe_internal,
                         FEValuesData<dim>& data) const;
     
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * FiniteElement.
-                                     */
     virtual void
     fill_fe_subface_values (const Mapping<dim> &mapping,
                            const typename Triangulation<dim>::cell_iterator &cell,
@@ -570,6 +542,112 @@ class FE_RaviartThomas : public FiniteElement<dim>
     template <int dim1> friend class FE_RaviartThomas;
 };
 
+
+
+/**
+ * The Raviart-Thomas elements with node functionals defined as point
+ * values in Gauss points.
+ *
+ * <h3>Description of node values</h3>
+ *
+ * For this Raviart-Thomas element, the node values are not cell and
+ * face moments with respect to certain polynomials, but the values in
+ * quadrature points.
+ *
+ * For an RT-element of degree <i>k</i>, we choose
+ * <i>k+1<sup>d-1</sup></i> Gauss points on each face. This way, the
+ * normal component which is in <i>Q<sub>k</sub></i> is uniquely
+ * determined. Furthermore, since this Gauss-formula is exact on
+ * <i>Q<sub>2k+1</sub></i>, these node values correspond to the exact
+ * integration of the moments of the RT-space.
+ *
+ * In the interior of the cells, the moments are with respect to an
+ * anisotropic <i>Q<sub>k</sub></i> space, where the test functions
+ * are one degree lower in the direction corresponding to the vector
+ * component under consideration. This can be emulated by using an
+ * anisotropic Gauss formula for integration.
+ *
+ * @warning The degree stored in the member variable
+ * FiniteElementData<dim>::degree is higher by one than the
+ * constructor argument!
+ * 
+ * @author Guido Kanschat, 2005
+ */
+template <int dim>
+class FE_RaviartThomasNodal
+  :
+  public FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim>
+{
+                                    /**
+                                     * Constructor for the Raviart-Thomas
+                                     * element of degree @p p.
+                                     */
+    FE_RaviartThomasNodal (const unsigned int p);
+    
+                                    /**
+                                     * Return a string that uniquely
+                                     * identifies a finite
+                                     * element. This class returns
+                                     * <tt>FE_RaviartThomasNodal<dim>(degree)</tt>, with
+                                     * @p dim and @p degree
+                                     * replaced by appropriate
+                                     * values.
+                                     */
+    virtual std::string get_name () const;
+
+    virtual FiniteElement<dim>* clone () const;
+
+                                    /**
+                                     * Check whether a shape function
+                                     * may be non-zero on a face.
+                                     *
+                                     * Right now, always returns
+                                     * @p true.
+                                     */
+    virtual bool has_support_on_face (const unsigned int shape_index,
+                                     const unsigned int face_index) const;    
+  private:
+                                    /**
+                                     * Only for internal use. Its
+                                     * full name is
+                                     * @p get_dofs_per_object_vector
+                                     * function and it creates the
+                                     * @p dofs_per_object vector that is
+                                     * needed within the constructor to
+                                     * be passed to the constructor of
+                                     * @p FiniteElementData.
+                                     */
+    static std::vector<unsigned int>
+    get_dpo_vector (const unsigned int degree);
+
+                                    /**
+                                     * Compute the vector used for
+                                     * the
+                                     * @p restriction_is_additive
+                                     * field passed to the base
+                                     * class's constructor.
+                                     */
+    static std::vector<bool>
+    get_ria_vector (const unsigned int degree);
+                                    /**
+                                     * Initialize the
+                                     * FiniteElementBase<dim>::unit_support_points
+                                     * and FiniteElementBase<dim>::unit_face_support_points
+                                     * fields. Called from the
+                                     * constructor.
+                                     */
+    void initialize_unit_support_points (const unsigned int degree);
+
+                                    /**
+                                     * Initialize the
+                                     * #inverse_node_matrix
+                                     * field. Called from the
+                                     * constructor.
+                                     */
+    void initialize_node_matrix ();
+};
+
+
 /*@}*/
 
 /* -------------- declaration of explicit specializations ------------- */
index 4c9cebdd16847e62deda3b898ea265f29359fab9..b6fb9f3a8f2a7f84e176b0afcf4ea077220eeb97 100644 (file)
@@ -26,7 +26,6 @@ FE_PolyTensor<POLY,dim>::FE_PolyTensor (
                FiniteElement<dim> (fe_data,
                                    restriction_is_additive_flags,
                                    nonzero_components),
-                degree(degree),
                 poly_space(POLY(degree))
 {
   cached_point(0) = -1;
@@ -252,12 +251,13 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
 
 template <class POLY, int dim>
 void
-FE_PolyTensor<POLY,dim>::fill_fe_values (const Mapping<dim>                   &mapping,
-                                  const typename Triangulation<dim>::cell_iterator &cell,
-                                  const Quadrature<dim>                &quadrature,
-                                  typename Mapping<dim>::InternalDataBase &mapping_data,
-                                  typename Mapping<dim>::InternalDataBase &fedata,
-                                  FEValuesData<dim>                    &data) const
+FE_PolyTensor<POLY,dim>::fill_fe_values (
+  const Mapping<dim>                   &mapping,
+  const typename Triangulation<dim>::cell_iterator &cell,
+  const Quadrature<dim>                &quadrature,
+  typename Mapping<dim>::InternalDataBase &mapping_data,
+  typename Mapping<dim>::InternalDataBase &fedata,
+  FEValuesData<dim>                    &data) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -269,7 +269,8 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (const Mapping<dim>                   &m
   const unsigned int n_quad = quadrature.n_quadrature_points;
   const UpdateFlags flags(fe_data.current_update_flags());
   
-  Assert(mapping_type == independent, ExcNotImplemented());
+  Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+        ExcNotImplemented());
   
   Assert(!(flags & update_values) || fe_data.shape_values.size() == n_dofs,
         ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs));
@@ -307,13 +308,14 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (const Mapping<dim>                   &m
 
 template <class POLY, int dim>
 void
-FE_PolyTensor<POLY,dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
-                               const typename Triangulation<dim>::cell_iterator &cell,
-                               const unsigned int                    face,
-                               const Quadrature<dim-1>              &quadrature,
-                               typename Mapping<dim>::InternalDataBase       &mapping_data,
-                               typename Mapping<dim>::InternalDataBase       &fedata,
-                               FEValuesData<dim>                    &data) const
+FE_PolyTensor<POLY,dim>::fill_fe_face_values (
+  const Mapping<dim>                   &mapping,
+  const typename Triangulation<dim>::cell_iterator &cell,
+  const unsigned int                    face,
+  const Quadrature<dim-1>              &quadrature,
+  typename Mapping<dim>::InternalDataBase       &mapping_data,
+  typename Mapping<dim>::InternalDataBase       &fedata,
+  FEValuesData<dim>                    &data) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -334,7 +336,8 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (const Mapping<dim>
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
-  Assert(mapping_type == independent, ExcNotImplemented());
+  Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+        ExcNotImplemented());
 //TODO: Size assertions
   
   for (unsigned int i=0; i<n_dofs; ++i)
@@ -365,14 +368,15 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (const Mapping<dim>
 
 template <class POLY, int dim>
 void
-FE_PolyTensor<POLY,dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
-                                  const typename Triangulation<dim>::cell_iterator &cell,
-                                  const unsigned int                    face,
-                                  const unsigned int                    subface,
-                                  const Quadrature<dim-1>              &quadrature,
-                                  typename Mapping<dim>::InternalDataBase       &mapping_data,
-                                  typename Mapping<dim>::InternalDataBase       &fedata,
-                                  FEValuesData<dim>                    &data) const
+FE_PolyTensor<POLY,dim>::fill_fe_subface_values (
+  const Mapping<dim>                   &mapping,
+  const typename Triangulation<dim>::cell_iterator &cell,
+  const unsigned int                    face,
+  const unsigned int                    subface,
+  const Quadrature<dim-1>              &quadrature,
+  typename Mapping<dim>::InternalDataBase       &mapping_data,
+  typename Mapping<dim>::InternalDataBase       &fedata,
+  FEValuesData<dim>                    &data) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -393,7 +397,8 @@ FE_PolyTensor<POLY,dim>::fill_fe_subface_values (const Mapping<dim>
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
-  Assert(mapping_type == independent, ExcNotImplemented());
+  Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+        ExcNotImplemented());
 //TODO: Size assertions
   
   for (unsigned int i=0; i<n_dofs; ++i)
@@ -450,5 +455,40 @@ FE_PolyTensor<POLY,dim>::element_multiplicity (const unsigned int index) const
 }
 
 
+template <class POLY, int dim>
+UpdateFlags
+FE_PolyTensor<POLY,dim>::update_once (const UpdateFlags flags) const
+{
+  Assert (mapping_type != no_mapping, ExcNotInitialized());
+  const bool values_once = (mapping_type == independent);
+  
+  UpdateFlags out = update_default;
+
+  if (values_once && (flags & update_values))
+    out |= update_values;
+
+  return out;
+}
+
+
+template <class POLY, int dim>
+UpdateFlags
+FE_PolyTensor<POLY,dim>::update_each (const UpdateFlags flags) const
+{
+  Assert (mapping_type != no_mapping, ExcNotInitialized());
+  const bool values_once = (mapping_type == independent);
+  
+  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_second_derivatives)
+    out |= update_second_derivatives | update_covariant_transformation;
+
+  return out;
+}
+
 
 template class FE_PolyTensor<PolynomialsRaviartThomas<deal_II_dimension>,deal_II_dimension>;
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
new file mode 100644 (file)
index 0000000..1570bf7
--- /dev/null
@@ -0,0 +1,294 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003, 2004, 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <base/quadrature_lib.h>
+#include <base/table.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/mapping.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_values.h>
+#include <fe/fe_tools.h>
+
+#ifdef HAVE_STD_STRINGSTREAM
+#  include <sstream>
+#else
+#  include <strstream>
+#endif
+
+
+template <int dim>
+FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
+               :
+               FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim> (
+                 deg,
+                 FiniteElementData<dim>(get_dpo_vector(deg),
+                                        dim, deg+1),
+                 get_ria_vector (deg),
+                 std::vector<std::vector<bool> >(
+                   FiniteElementData<dim>(get_dpo_vector(deg),
+                                          dim,deg+1).dofs_per_cell,
+                   std::vector<bool>(dim,true)))
+{
+  Assert (dim >= 2, ExcImpossibleInDim(dim));
+  
+  this->mapping_type = independent_on_cartesian;
+  
+  for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
+    this->prolongation[i].reinit (this->dofs_per_cell,
+                                 this->dofs_per_cell);
+  FETools::compute_embedding_matrices (*this, &this->prolongation[0]);
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points (deg);
+  initialize_node_matrix();
+}
+
+
+
+template <int dim>
+std::string
+FE_RaviartThomasNodal<dim>::get_name () const
+{
+                                  // note that the
+                                  // FETools::get_fe_from_name
+                                  // function depends on the
+                                  // particular format of the string
+                                  // this function returns, so they
+                                  // have to be kept in synch
+
+#ifdef HAVE_STD_STRINGSTREAM
+  std::ostringstream namebuf;
+#else
+  std::ostrstream namebuf;
+#endif
+  
+  namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << degree-1 << ")";
+
+#ifndef HAVE_STD_STRINGSTREAM
+  namebuf << std::ends;
+#endif
+  return namebuf.str();
+}
+
+
+
+template <int dim>
+bool
+FE_RaviartThomasNodal<dim>::has_support_on_face (unsigned int, unsigned int) const
+{
+  return true;
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_RaviartThomasNodal<dim>::clone() const
+{
+  return new FE_RaviartThomasNodal<dim>(degree-1);
+}
+
+
+#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=0; d<dim-1; ++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)
+{
+  unsigned int dofs_per_cell, dofs_per_face;
+  switch (dim)
+    {
+      case 2:
+           dofs_per_face = deg+1;
+           dofs_per_cell = 2*(deg+1)*(deg+2);
+           break;
+      case 3:
+           dofs_per_face = (deg+1)*(deg+1);
+           dofs_per_cell = 3*(deg+1)*(deg+1)*(deg+2);
+           break;
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+  Assert (FiniteElementData<dim>(get_dpo_vector(deg),dim).dofs_per_cell ==
+         dofs_per_cell,
+         ExcInternalError());
+  Assert (FiniteElementData<dim>(get_dpo_vector(deg),dim).dofs_per_face ==
+         dofs_per_face,
+         ExcInternalError());
+  
+                                  // 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_unit_support_points (const unsigned int deg)
+{
+  this->unit_support_points.resize (this->dofs_per_cell);
+  this->unit_face_support_points.resize (this->dofs_per_face);
+  
+  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.n_quadrature_points == this->dofs_per_face,
+             ExcInternalError());
+      for (unsigned int k=0;k<this->dofs_per_face;++k)
+       this->unit_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<faces.n_quadrature_points;++k)
+       this->unit_support_points[k] = faces.point(k);
+
+      current = faces.n_quadrature_points;
+    }
+                                  // 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));
+      for (unsigned int k=0;k<quadrature->n_quadrature_points;++k)
+       this->unit_support_points[current++] = quadrature->point(k);
+      delete quadrature;
+    }
+  Assert (current == this->dofs_per_cell, ExcInternalError());
+}
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::initialize_node_matrix ()
+{
+  const unsigned int n_dofs = this->dofs_per_cell;
+
+                                  // 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> N(n_dofs, n_dofs);
+
+                                  // The curent node functional index
+  unsigned int current = 0;
+
+                                  // For each face and all quadrature
+                                  // points on it, the node value is
+                                  // the normal component of the
+                                  // shape function, possibly
+                                  // pointing in negative direction.
+  for (unsigned int face = 0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int k=0;k<this->dofs_per_face;++k)
+      {
+       for (unsigned int i=0;i<n_dofs;++i)
+         N(current,i) = this->shape_value_component(
+           i, unit_support_point(current),
+           GeometryInfo< dim >::unit_normal_direction[face])
+                        * GeometryInfo< dim >::unit_normal_orientation[face];
+       ++current;
+      }
+                                  // Interior degrees of freedom in each direction
+  const unsigned int n_cell = (n_dofs - current) / dim;
+  
+  for (unsigned int d=0;d<dim;++d)
+    for (unsigned int k=0;k<n_cell;++k)
+      {
+       for (unsigned int i=0;i<n_dofs;++i)
+         N(current,i) = this->shape_value_component(i, unit_support_point(current), d);
+       ++current;
+      }
+  Assert (current == n_dofs, ExcInternalError());
+
+  inverse_node_matrix.invert(N);
+}
+
+
+
+template FE_RaviartThomasNodal<deal_II_dimension>;

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.