]> https://gitweb.dealii.org/ - dealii.git/commitdiff
BDM elements
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 18 Jul 2010 05:02:58 +0000 (05:02 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 18 Jul 2010 05:02:58 +0000 (05:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@21515 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/deal.II/include/fe/fe_bdm.h b/deal.II/deal.II/include/fe/fe_bdm.h
new file mode 100644 (file)
index 0000000..228c09f
--- /dev/null
@@ -0,0 +1,124 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 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_bdm_h
+#define __deal2__fe_bdm_h
+
+#include <base/config.h>
+#include <base/table.h>
+#include <base/polynomials_bdm.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>
+
+using namespace dealii;
+
+/**
+ * The Brezzi-Douglas-Marini element.
+ *
+ * <h3>Degrees of freedom</h3>
+ *
+ * @todo This is for 2D only.
+ *
+ * @todo Transformation works only for uniform, Cartesian meshes
+ *
+ * The BDM element of order @p p has <i>p+1</i> degrees of freedom on
+ * each face. These are implemented as the function values in the
+ * <i>p+1</i> Gauss points on each face.
+ *
+ * Additionally, for order greater or equal 2, we have additional
+ * <i>p(p-1)<i>, the number of vector valued polynomials in
+ * <i>P<sub>p</sub></i>, interior degrees of freedom. These are the
+ * vector function values in the first <i>p(p-1)/2<i> of the
+ * <i>p<sup>2</sup></i> Gauss points in the cell.
+ */
+template <int dim>
+class FE_BDM
+  :
+  public FE_PolyTensor<PolynomialsBDM<dim>, dim>
+{
+  public:
+                                    /**
+                                     * Constructor for the BDM
+                                     * element of degree @p p.
+                                     */
+    FE_BDM (const unsigned int p);
+    
+                                    /**
+                                     * Return a string that uniquely
+                                     * identifies a finite
+                                     * element. This class returns
+                                     * <tt>FE_BDM<dim>(degree)</tt>, with
+                                     * @p dim and @p degree
+                                     * replaced by appropriate
+                                     * values.
+                                     */
+    virtual std::string get_name () const;
+    
+    virtual FiniteElement<dim>* clone () 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;
+  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
+                                     * FiniteElement<dim>::generalized_support_points
+                                     * and FiniteElement<dim>::generalized_face_support_points
+                                     * fields. Called from the
+                                     * constructor.
+                                     */
+    void initialize_support_points (const unsigned int rt_degree);
+                                    /**
+                                     * The values in the interior
+                                     * support points of the
+                                     * polynomials needed as test
+                                     * functions. The outer vector is
+                                     * indexed by quadrature points,
+                                     * the inner by the test
+                                     * function.
+                                     */
+    std::vector<std::vector<double> > test_values;
+};
+
+#endif
diff --git a/deal.II/deal.II/source/fe/fe_bdm.cc b/deal.II/deal.II/source/fe/fe_bdm.cc
new file mode 100644 (file)
index 0000000..1822c27
--- /dev/null
@@ -0,0 +1,375 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    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
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <base/quadrature_lib.h>
+#include <base/qprojector.h>
+#include <base/polynomials_p.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_bdm.h>
+#include <fe/fe_values.h>
+#include <fe/fe_tools.h>
+
+#include <iostream>
+#include <sstream>
+
+template <int dim>
+FE_BDM<dim>::FE_BDM (const unsigned int deg)
+               :
+               FE_PolyTensor<PolynomialsBDM<dim>, dim> (
+                 deg,
+                 FiniteElementData<dim>(get_dpo_vector(deg),
+                                        dim, deg+1, FiniteElementData<dim>::Hdiv, 1),
+                 get_ria_vector (deg),
+                 std::vector<std::vector<bool> >(PolynomialsBDM<dim>::compute_n_pols(deg),
+                                                 std::vector<bool>(dim,true)))
+{
+  Assert (dim >= 2, ExcImpossibleInDim(dim));
+  const unsigned int n_dofs = this->dofs_per_cell;
+  
+  this->mapping_type = mapping_piola;
+                                  // 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
+                                  //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);
+
+//   std::cout << std::endl;
+//   M.print_formatted(std::cout, 2, true);
+  
+  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
+  this->inverse_node_matrix.invert(M);
+                                  // From now on, the shape functions
+                                  // will be the correct ones, not
+                                  // the raw shape functions anymore.
+  
+  this->reinit_restriction_and_prolongation_matrices(true, true);
+  FETools::compute_embedding_matrices (*this, this->prolongation, true);
+  
+//   FullMatrix<double> face_embeddings[GeometryInfo<dim>::subfaces_per_face];
+//   for (unsigned int i=0; i<GeometryInfo<dim>::subfaces_per_face; ++i)
+//     face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
+//   FETools::compute_face_embedding_matrices(*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;
+//   for (unsigned int d=0;d<GeometryInfo<dim>::subfaces_per_face;++d)
+//     for (unsigned int i=0;i<face_embeddings[d].m();++i)
+//       {
+//     for (unsigned int j=0;j<face_embeddings[d].n();++j)
+//       this->interface_constraints(target_row,j) = face_embeddings[d](i,j);
+//     ++target_row;
+//       }
+}
+
+
+
+template <int dim>
+std::string
+FE_BDM<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
+
+                                  // 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_BDM<" << dim << ">(" << this->degree-1 << ")";
+
+  return namebuf.str();
+}
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_BDM<dim>::clone() const
+{
+  return new FE_BDM<dim>(*this);
+}
+
+
+
+template <int dim>
+void
+FE_BDM<dim>::interpolate(
+  std::vector<double>&,
+  const std::vector<double>&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <int dim>
+void
+FE_BDM<dim>::interpolate(
+  std::vector<double>&,
+  const std::vector<Vector<double> >&,
+  unsigned int) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+
+template <int dim>
+void
+FE_BDM<dim>::interpolate(
+  std::vector<double>& local_dofs,
+  const VectorSlice<const std::vector<std::vector<double> > >& values) const
+{
+  AssertDimension (values.size(), dim);
+  Assert (values[0].size() == this->generalized_support_points.size(),
+         ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
+  Assert (local_dofs.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+                                  // First do interpolation on
+                                  // faces. There, the component
+                                  // evaluated depends on the face
+                                  // direction and orientation.
+  unsigned int fbase = 0;
+  unsigned int f=0;
+  for (;f<GeometryInfo<dim>::faces_per_cell;
+       ++f, fbase+=this->dofs_per_face)
+    {
+      for (unsigned int i=0;i<this->dofs_per_face;++i)
+       {
+         local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i];
+       }
+    }
+
+                                  // Done for BDM1
+  if (fbase == this->dofs_per_cell) return;
+  
+                                  // What's missing are the interior
+                                  // degrees of freedom. In each
+                                  // point, we take all components of
+                                  // the solution.
+  Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+                                  // Here, the number of the point
+                                  // and of the shape function
+                                  // coincides. This will change
+                                  // below, since we have more
+                                  // support points than test
+                                  // functions in the interior.
+  const unsigned int pbase = fbase;
+  for (unsigned int d=0;d<dim;++d, fbase += test_values[0].size())
+    {
+      for (unsigned int i=0;i<test_values[0].size();++i)
+       {
+         local_dofs[fbase+i] = 0.;
+         for (unsigned int k=0;k<test_values.size();++k)
+           local_dofs[fbase+i] += values[d][pbase+k] * test_values[k][i];
+       }
+    }
+    
+  Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<unsigned int>
+FE_BDM<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_BDM<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
+  unsigned int
+    interior_dofs = dim*deg*(deg-1)/2;
+  if (dim>2)
+    {
+      interior_dofs *= deg-2;
+      interior_dofs /= 3;
+    }
+  
+  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_BDM<1>::get_ria_vector (const unsigned int)
+{
+  Assert (false, ExcImpossibleInDim(1));
+  return std::vector<bool>();
+}
+
+#endif
+
+
+template <int dim>
+std::vector<bool>
+FE_BDM<dim>::get_ria_vector (const unsigned int deg)
+{
+  Assert(dim==2, ExcNotImplemented());
+  const unsigned int dofs_per_cell = PolynomialsBDM<dim>::compute_n_pols(deg);
+  const unsigned int dofs_per_face = deg+1;
+                                  // all 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_BDM<dim>::initialize_support_points (const unsigned int deg)
+{
+                                  // interior point in 1d
+  unsigned int npoints = deg;
+                                  // interior point in 2d
+  if (dim >= 2)
+    {
+      npoints *= deg;
+//      npoints /= 2;
+    }
+                                  // interior point in 2d
+  if (dim >= 3)
+    {
+      npoints *= deg;
+//      npoints /= 3;
+    }
+  npoints += GeometryInfo<dim>::faces_per_cell * this->dofs_per_face;
+  
+  this->generalized_support_points.resize (npoints);
+  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 BDM 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<=1) return;
+                                  // Although the polynomial space is
+                                  // only P_{k-2}, we use the tensor
+                                  // product points for Q_{k-2}
+  QGauss<dim> quadrature(deg);
+  
+                                  // Remember where interior points start
+  const unsigned int ibase=current;
+//  for (unsigned int k=0;k<deg-1;++k)
+  for (unsigned int j=0;j<deg;++j)
+    for (unsigned int i=0;i<deg;++i)
+      {
+       this->generalized_support_points[current] = quadrature.point(current-ibase);
+       ++current;
+      }
+  Assert(current == npoints, ExcInternalError());
+
+                                  // Finaly, compute the values of
+                                  // the test functios in the
+                                  // interior quadrature points
+  PolynomialsP<dim> poly(deg-2);
+  
+  test_values.resize(quadrature.size());
+  std::vector<Tensor<1,dim> > dummy1;
+  std::vector<Tensor<2,dim> > dummy2;
+  
+  for (unsigned int k=0;k<quadrature.size();++k)
+    {
+      test_values[k].resize(poly.n());
+      poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2);
+      for (unsigned int i=0; i < poly.n(); ++i)
+       {
+         test_values[k][i] *= quadrature.weight(k);
+       }
+    }
+}
+
+
+template class FE_BDM<deal_II_dimension>;
+
index c7d3e62f95dd85edebb7f8dcf72173a167ddc57a..98c4ba1618afcf20dac1daa2eb5bed5ae90227f9 100644 (file)
@@ -150,6 +150,12 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+
+  <li><p> New: Brezzi-Douglas-Marini elements of arbitrary order in FE_BDM.
+  <br>
+  (GK 2010/07/19)
+  </p>
+  
   <li>
   <p>
   Fixed: The FEValues::get_cell() function was unusable from user code

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.