]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new vector valued DG elements
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 17 Sep 2010 17:07:25 +0000 (17:07 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 17 Sep 2010 17:07:25 +0000 (17:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@22018 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_dg_vector.h [new file with mode: 0644]
deal.II/deal.II/include/fe/fe_dg_vector.templates.h [new file with mode: 0644]
deal.II/deal.II/include/fe/fe_tools.h
deal.II/deal.II/source/fe/fe_dg_vector.cc [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/doc/news/changes.h

diff --git a/deal.II/deal.II/include/fe/fe_dg_vector.h b/deal.II/deal.II/include/fe/fe_dg_vector.h
new file mode 100644 (file)
index 0000000..0400d52
--- /dev/null
@@ -0,0 +1,186 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_raviart_thomas.h 18166 2009-01-09 21:21:16Z kanschat $
+//
+//    Copyright (C) 2010 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__fe_dg_vector_h
+#define __deal2__fe_dg_vector_h
+
+#include <base/config.h>
+#include <base/table.h>
+#include <base/polynomials_raviart_thomas.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <base/geometry_info.h>
+#include <fe/fe.h>
+#include <fe/fe_poly_tensor.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int spacedim> class MappingQ;
+
+
+/**
+ * DG elements based on vector valued polynomials.
+ *
+ * @ingroup fe
+ * @author Guido Kanschat
+ * @date 2010
+ */
+template <class POLY, int dim, int spacedim=dim>
+class FE_DGVector
+  :
+  public FE_PolyTensor<POLY, dim, spacedim>
+{
+  public:
+                                    /**
+                                     * Constructor for the vector
+                                     * element of degree @p p.
+                                     */
+    FE_DGVector (const unsigned int p, MappingType m);
+  public:
+    
+    FiniteElement<dim, spacedim>* clone() const;
+    
+                                    /**
+                                     * Return a string that uniquely
+                                     * identifies a finite
+                                     * element. This class returns
+                                     * <tt>FE_RaviartThomas<dim>(degree)</tt>, with
+                                     * @p dim and @p degree
+                                     * replaced by appropriate
+                                     * values.
+                                     */
+    virtual std::string get_name () const;
+
+
+                                    /**
+                                     * Check whether a shape function
+                                     * may be non-zero on a face.
+                                     *
+                                     * Returns always
+                                     * @p true.
+                                     */
+    virtual bool has_support_on_face (const unsigned int shape_index,
+                                     const unsigned int face_index) const;
+    
+    virtual void interpolate(std::vector<double>&                local_dofs,
+                            const std::vector<double>& values) const;
+    virtual void interpolate(std::vector<double>&                local_dofs,
+                            const std::vector<Vector<double> >& values,
+                            unsigned int offset = 0) const;
+    virtual void interpolate(
+      std::vector<double>& local_dofs,
+      const VectorSlice<const std::vector<std::vector<double> > >& values) const;
+    virtual unsigned int memory_consumption () 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);
+
+                                    /**
+                                     * Initialize the @p
+                                     * generalized_support_points
+                                     * field of the FiniteElement
+                                     * class and fill the tables with
+                                     * #interior_weights. Called
+                                     * from the constructor.
+                                     */
+    void initialize_support_points (const unsigned int degree);
+
+                                    /**
+                                     * Initialize the interpolation
+                                     * from functions on refined mesh
+                                     * cells onto the father
+                                     * cell. According to the
+                                     * philosophy of the
+                                     * Raviart-Thomas element, this
+                                     * restriction operator preserves
+                                     * the divergence of a function
+                                     * weakly.
+                                     */
+    void initialize_restriction ();
+    
+                                    /**
+                                     * Fields of cell-independent data.
+                                     *
+                                     * For information about the
+                                     * general purpose of this class,
+                                     * see the documentation of the
+                                     * base class.
+                                     */
+    class InternalData : public FiniteElement<dim>::InternalDataBase
+    {
+      public:
+                                        /**
+                                         * Array with shape function
+                                         * values in quadrature
+                                         * points. There is one row
+                                         * for each shape function,
+                                         * containing values for each
+                                         * quadrature point. Since
+                                         * the shape functions are
+                                         * vector-valued (with as
+                                         * many components as there
+                                         * are space dimensions), the
+                                         * value is a tensor.
+                                         *
+                                         * In this array, we store
+                                         * the values of the shape
+                                         * function in the quadrature
+                                         * points on the unit
+                                         * cell. The transformation
+                                         * to the real space cell is
+                                         * then simply done by
+                                         * multiplication with the
+                                         * Jacobian of the mapping.
+                                         */
+       std::vector<std::vector<Tensor<1,dim> > > shape_values;
+
+                                        /**
+                                         * Array with shape function
+                                         * gradients in quadrature
+                                         * points. There is one
+                                         * row for each shape
+                                         * function, containing
+                                         * values for each quadrature
+                                         * point.
+                                         *
+                                         * We store the gradients in
+                                         * the quadrature points on
+                                         * the unit cell. We then
+                                         * only have to apply the
+                                         * transformation (which is a
+                                         * matrix-vector
+                                         * multiplication) when
+                                         * visiting an actual cell.
+                                         */
+       std::vector<std::vector<Tensor<2,dim> > > shape_gradients;
+    };
+    Table<3, double> interior_weights;    
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/fe/fe_dg_vector.templates.h b/deal.II/deal.II/include/fe/fe_dg_vector.templates.h
new file mode 100644 (file)
index 0000000..46441f3
--- /dev/null
@@ -0,0 +1,106 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_poly.templates.h 21954 2010-09-13 20:49:06Z kanschat $
+//
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 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 <fe/fe_dg_vector.h>
+#include <fe/fe_tools.h>
+#include <base/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+//TODO:[GK] deg+1 is wrong here and should bew fixed after FiniteElementData was cleaned up
+
+template <class POLY, int dim, int spacedim>
+FE_DGVector<POLY,dim,spacedim>::FE_DGVector (
+  const unsigned int deg, MappingType map)
+               :
+               FE_PolyTensor<POLY, dim, spacedim>(
+                 deg,
+                 FiniteElementData<dim>(
+                   get_dpo_vector(deg), dim, deg+1, FiniteElementData<dim>::L2, 1),
+                 std::vector<bool>(POLY::compute_n_pols(deg), true),
+                 std::vector<std::vector<bool> >(POLY::compute_n_pols(deg),
+                                                 std::vector<bool>(dim,true)))
+{
+  this->mapping_type = map;
+  const unsigned int polynomial_degree = this->tensor_degree();
+  
+  QGauss<dim> quadrature(polynomial_degree+1);
+  this->generalized_support_points = quadrature.get_points();
+  
+ this->reinit_restriction_and_prolongation_matrices(true, true);
+ FETools::compute_projection_matrices (*this, this->restriction, true);
+ FETools::compute_embedding_matrices (*this, this->prolongation, true);
+}
+
+               
+template <class POLY, int dim, int spacedim>
+FiniteElement<dim, spacedim> *
+FE_DGVector<POLY,dim,spacedim>::clone() const
+{
+  return new FE_DGVector<POLY, dim, spacedim>(*this);
+}
+
+
+template <class POLY, int dim, int spacedim>
+std::string
+FE_DGVector<POLY,dim,spacedim>::get_name() const
+{
+  std::ostringstream namebuf;  
+  namebuf << "FE_DGVector_" << this->poly_space.name()
+         << "<" << dim << ">(" << this->degree-1 << ")";
+  return namebuf.str();
+}
+
+
+template <class POLY, int dim, int spacedim>
+std::vector<unsigned int>
+FE_DGVector<POLY,dim,spacedim>::get_dpo_vector (const unsigned int deg)
+{
+  std::vector<unsigned int> dpo(dim+1);
+  dpo[dim] = POLY::compute_n_pols(deg);
+  
+  return dpo;
+}
+
+
+template <class POLY, int dim, int spacedim>
+bool
+FE_DGVector<POLY,dim,spacedim>::has_support_on_face (
+  const unsigned int,
+  const unsigned int) const
+{
+  return true;
+}
+
+
+template <class POLY, int dim, int spacedim>
+void
+FE_DGVector<POLY,dim,spacedim>::interpolate(
+  std::vector<double>&,
+  const std::vector<double>&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <class POLY, int dim, int spacedim>
+void
+FE_DGVector<POLY,dim,spacedim>::interpolate(
+  std::vector<double>& /*local_dofs*/,
+  const std::vector<Vector<double> >& /*values*/,
+  unsigned int /*offset*/) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+DEAL_II_NAMESPACE_CLOSE
index 10860141a45bd1ba4e4919d5521ad69b4d3fb59b..acb9b373f917e342ce7707051d170a562976a331 100644 (file)
@@ -477,8 +477,10 @@ class FETools
                                      * to use this function mostly.
                                      */
     template <int dim, typename number, int spacedim>
-    static void compute_projection_matrices(const FiniteElement<dim,spacedim> &fe,
-                                           std::vector<std::vector<FullMatrix<number> > >& matrices);
+    static void compute_projection_matrices(
+      const FiniteElement<dim,spacedim> &fe,
+      std::vector<std::vector<FullMatrix<number> > >& matrices,
+      const bool isotropic_only = false);
 
                                     /**
                                       * Projects scalar data defined in
diff --git a/deal.II/deal.II/source/fe/fe_dg_vector.cc b/deal.II/deal.II/source/fe/fe_dg_vector.cc
new file mode 100644 (file)
index 0000000..dac4f7c
--- /dev/null
@@ -0,0 +1,27 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_dgp.cc 17866 2008-12-05 22:27:44Z bangerth $
+//
+//    Copyright (C) 2010 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 <fe/fe_dg_vector.templates.h>
+#include <base/polynomials_abf.h>
+#include <base/polynomials_bdm.h>
+#include <base/polynomials_nedelec.h>
+#include <base/polynomials_raviart_thomas.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template class FE_DGVector<PolynomialsABF<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsBDM<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsNedelec<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsRaviartThomas<deal_II_dimension>, deal_II_dimension>;
+
+DEAL_II_NAMESPACE_CLOSE
+
index 58cf5477412e275b35c404065fe65d743f75efeb..d9a6f3dfdf1304e8bb904f75c1f37dd43442423d 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -457,8 +457,9 @@ FE_RaviartThomas<dim>::get_dpo_vector (const unsigned int deg)
 
 template <int dim>
 bool
-FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
-                                            const unsigned int face_index) const
+FE_RaviartThomas<dim>::has_support_on_face (
+  const unsigned int shape_index,
+  const unsigned int face_index) const
 {
   Assert (shape_index < this->dofs_per_cell,
          ExcIndexRange (shape_index, 0, this->dofs_per_cell));
index 7a9d4320c1ac70916500d8e5d4876c0e03c4fadc..76082357d00f7dbbc89d146771099296d4d3ae83 100644 (file)
@@ -933,7 +933,7 @@ FETools::compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
 template <>
 void
 FETools::compute_projection_matrices(const FiniteElement<1,2>&,
-                                    std::vector<std::vector<FullMatrix<double> > >& )
+                                    std::vector<std::vector<FullMatrix<double> > >&, bool)
 {
   Assert(false, ExcNotImplemented());
 }
@@ -942,7 +942,7 @@ FETools::compute_projection_matrices(const FiniteElement<1,2>&,
 template <>
 void
 FETools::compute_projection_matrices(const FiniteElement<2,3>&,
-                                    std::vector<std::vector<FullMatrix<double> > >& )
+                                    std::vector<std::vector<FullMatrix<double> > >&, bool)
 {
   Assert(false, ExcNotImplemented());
 }
@@ -953,7 +953,8 @@ FETools::compute_projection_matrices(const FiniteElement<2,3>&,
 template <int dim, typename number, int spacedim>
 void
 FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
-                                    std::vector<std::vector<FullMatrix<number> > >& matrices)
+                                    std::vector<std::vector<FullMatrix<number> > >& matrices,
+                                    const bool isotropic_only)
 {
   const unsigned int n  = fe.dofs_per_cell;
   const unsigned int nd = fe.n_components();
@@ -1005,9 +1006,12 @@ FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
     mass.gauss_jordan();
   }
 
-                                  // loop over all possible refinement cases
-  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
-       ref_case < RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
+                                  // loop over all possible
+                                  // refinement cases
+  unsigned int ref_case = (isotropic_only)
+                         ? RefinementCase<dim>::isotropic_refinement
+                         : RefinementCase<dim>::cut_x;
+  for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
     {
       const unsigned int
        nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
@@ -2257,7 +2261,7 @@ void FETools::compute_face_embedding_matrices<deal_II_dimension,double>
 
 template
 void FETools::compute_projection_matrices<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&);
+(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&, bool);
 
 template
 void FETools::interpolate<deal_II_dimension>
index 0b75000b1ecd023d70936255b3df181030bffaaf..ceadb44286638836d8c877ec6288bb00d9d00f87 100644 (file)
@@ -267,6 +267,13 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
 <h3>deal.II</h3>
 
 <ol>
+
+  <li><p>New: FE_DGVector implements discontinuous elements based on
+  vector valued polynomial spaces.
+  <br>
+  (GK 2010/09/17)
+  </p></li>
+
   <li>
   <p>
   Fixed: The methods VectorTools::interpolate_boundary_values and

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.